bugGNU Octave - Bugs: bug #44004, eig and eigs routine

 
 

bug #44004: eig and eigs routine

Submitter:  Sebastian Glane <sebastian_g>
Submitted:  Wed 14 Jan 2015 06:06:42 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Inaccurate Result
Status:  Invalid / Not an Octave Bug Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * 3.8.1 Operating System:  * GNU/Linux
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Wed 25 Nov 2015 12:58:31 PM UTC, comment #9: 

The original comment does not specified which version of ARPACK is involved. If 3.1.5, dneupd is buggy.

Marco

Marco Caliari <caliari>
Group Member
Wed 25 Nov 2015 12:19:44 AM UTC, comment #8: 

Can this bug be closed as not an Octave bug, but an upstream bug in ARPACK's implementation of dneupd?

Rik <rik5>
Group administrator
Fri 13 Feb 2015 11:03:52 AM UTC, comment #7: 

Dear all,

there is a bug in ARPACK up to 3.1.5 related to dneupd. It was fixed in 3.2.0 and I'm proposing a better fix here https://github.com/opencollab/arpack-ng/pull/2

It would be nice if you can check that it really fix your problems.

Marco

Marco Caliari <caliari>
Group Member
Thu 29 Jan 2015 05:30:38 AM UTC, comment #6: 

I have just thought of an ugly work-around for this problem.  Before I create the patch, could someone (Jordi?) comment on whether or not it would be suitable?

It seems that  dneupd  returns the eigenvalues in decreasing order when it works correctly.  However, when it fails the results are not monotonic (since they are the largest eigenvalues).  We could check if the results are monotonic, and if not, repeat the calculation.

Yes, I said it is ugly...  Thoughts?

Lachlan Andrew <lachlan>
Sat 24 Jan 2015 05:41:13 AM UTC, comment #5: 

Thanks, Jordi, for pointing out that the initialization is random.  However, that is not causing the problem here.

It seems that the call to the   dnaupd   function in ARPACK works reliably (and that is where the randomness in the algorithm is).

The problem is in the call to  dneupd.  Most of the time, it correctly returns the first  k  values of workl[iptr(5)] as expected.  However sometimes it returns the last two (which are not valid eigenvalues) instead of the first two.

Sometimes it replaces three or four, depending on the value of  k.

Note that the main difference between
  D = eigs(...)
and
  [V,D] = eigs(...)
is in the first parameter passed to  dneupd.

bug #41305 also suggests that there is a (different) bug in  dneupd, so I wouldn't rule out that possibility.

Lachlan Andrew <lachlan>
Wed 21 Jan 2015 02:33:19 PM UTC, comment #4: 

As you can read  [here](http://scicomp.stackexchange.com/questions/18688/eigs-routine-in-octave) , see the comments, the error does not occur if eigs is called with

D=eigs(...)

instead of

[V,D]=eigs(...).

Therefore, to me it seems like an octave bug not a problem with ARPACK. But I am not developer...

Sebastian Glane <sebastian_g>
Tue 20 Jan 2015 09:09:28 PM UTC, comment #3: 


> The function eigs should not return different results each time it is run.


The Arnoldi method involves some initial randomness. Getting slightly different results each time is normal:

https://en.wikipedia.org/wiki/Arnoldi_iteration#Krylov_subspaces_and_the_power_iteration

Jordi GutiƩrrez Hermoso <jordigh>
Group Member
Tue 20 Jan 2015 06:14:53 PM UTC, comment #2: 

I modified the code to ensure convergence (flag -- fl) and to check if the eigenvalues are real. If eigenvalues may be complex the sort command can cause the incorrect results.

% solves the Eigenvalue Problem
%   -y''= lambda^2 y on the interval (a,b)
%
% with mixed boundary conditions
%   -- dirichlet boundary condition y(a)=0
%   -- robin boundary condition   y'(b)=y(b)
clear all;
close all;
p=6;
%% number of points
n=2^p;
%% grid space
h=1/2^p;
%% numbers of eigenvalues to find
kmax=8;
%% assemble symmetric 2nd finite difference matrix
Lh=-gallery('tridiag',n,1,-2,1)./(h^2);
%% apply robin boundary in last line
Lh(n,:)=0.0;
Lh(n,n-1)=-2/h^2;
Lh(n,n)=-1/h^2*(2*h-2);
%% options for eigs
% opts.tol=1e-10;
opts.p=max(floor(n/2),2*kmax);
comp=0;
while comp<10 % do the same thing 10 times
%% use eig to get all eigenvalues
[~,LambdaEig]=eig(Lh);
if isreal(diag(LambdaEig))
LambdaEig=sort(diag(LambdaEig));
lambdaEig=LambdaEig(1:kmax);
else
warning('eigenvalues of eig are not real')
end
%% use eigs routine to get 1st kmax eigenvalues of smallest magnitude
[~,LambdaEigs,fl]=eigs(Lh,kmax,'sr',opts);
if isreal(diag(LambdaEigs))
LambdaEig=sort(diag(LambdaEig));
lambdaEigs=sort(diag(LambdaEigs));
else
warning('eigenvalues of eigs are not real')
end
%% compute difference of 1st kmax eigenvalues of smallest magnitude
errmax=max(abs(lambdaEig-lambdaEigs));
disp(['eigs flag: ' num2str(fl) ', maximum error of eigenvalues: '...
num2str(errmax)])
comp=comp+1;
end

Sebastian Glane <sebastian_g>
Tue 20 Jan 2015 02:15:15 AM UTC, comment #1: 

The function  eigs  should not return different results each time it is run.  It is clearly returning incorrect results.  However, this function uses a third-party library; we'll have to see if the bug is in that library or Octave's own code.

Lachlan Andrew <lachlan>
Wed 14 Jan 2015 06:06:42 PM UTC, original submission:  

The code below demonstrates a large difference in the computed eigenvalues of the eig- and eigs-routine. eigs will find complex pairs whereas eig finds real ones.

It gets much better if you replace 'sm' by 'sr'. Execute the script a few times and observe that somehow large difference occure. This seems unreliable to me.

By the way, if you choose 'sr' and p=3, you will get the error

sigma = "sr" requires complex or unsymmetric problem



% solves Eigenvalue Problem
%   -y''= lambda^2 y on the interval (a,b)
%
% with mixed boundary conditions
%   -- dirichlet boundary condition y(a)=0
%   -- robin boundary condition   y'(b)=y(b)
clear all;
close all;
p=6;
%% number of points
n=2^p;
%% grid space
h=1/2^p;
%% numbers of eigenvalues to find
kmax=8;
%% assemble symmetric 2nd finite difference matrix
Lh=-gallery('tridiag',n,1,-2,1)./(h^2);
%% apply robin boundary in last line
Lh(n,:)=0.0;
Lh(n,n-1)=-2/h^2;
Lh(n,n)=-1/h^2*(2*h-2);
%% use eig to get all eigenvalues
[~,LambdaEig]=eig(Lh);
LambdaEig=sort(diag(LambdaEig));
lambdaEig=LambdaEig(1:kmax);
%% use eigs routine to get 1st kmax eigenvalues of smallest magnitude
[~,LambdaEigs,fl]=eigs(Lh,kmax,'sm');
lambdaEigs=sort(diag(LambdaEigs))(1:kmax);
%% compute difference of 1st kmax eigenvalues of smallest magnitude
err=abs(lambdaEig-lambdaEigs)


Sebastian Glane <sebastian_g>

 

(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

Digest:
   bug dependencies.

 

Carbon-Copy List
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by caliari (Posted a comment)
  • -email is unavailable- added by jordigh (Posted a comment)
  • -email is unavailable- added by lachlan (Posted a comment)
  • -email is unavailable- added by sebastian_g (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 4 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2016-02-22 mtmiller Dependencies- bugs #45153 is dependent
    2015-11-25 rik5 StatusNeed Info Invalid / Not an Octave Bug
        Open/ClosedOpen Closed
    2015-11-25 rik5 StatusNone Need Info

    Back to the top

    Powered by Savane 3.13-f8d8.
    Corresponding source code