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

 
 

bug #44004: eig and eigs routine

Submitted by:  Sebastian Glane <sebastian_g>
Submitted on:  Wed 14 Jan 2015 06:06:42 PM UTC  
 
Category: Octave FunctionSeverity: 3 - Normal
Priority: 5 - NormalItem Group: Inaccurate Result
Status: InvalidAssigned to: None
Originator Name: Open/Closed: Closed
Release: 3.8.1Operating System: GNU/Linux

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 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>
Project 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>
Project 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>
Project 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>
Project Member
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>
Project Member
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>
Project Administrator
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(2h-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>
Project Member
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

Sebastian Glane <sebastian_g>

 

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

Attach File(s):
   
   
Comment:
   

No files currently attached

 

Depends on the following items: None found

Digest:
   bug dependencies.

 

Carbon-Copy List
  • -unavailable- added by rik5 (Posted a comment)
  • -unavailable- added by caliari (Posted a comment)
  • -unavailable- added by jordigh (Posted a comment)
  • -unavailable- added by lachlan (Posted a comment)
  • -unavailable- added by sebastian_g (Submitted the item)
  •  

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

    Only project members can vote.

     

    Please enter the title of George Orwell's famous dystopian book (it's a date):

     

     

    Follow 4 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Mon 22 Feb 2016 12:44:02 AM UTCmtmillerDependencies-=>bugs #45153 is dependent
    Wed 25 Nov 2015 04:13:10 PM UTCrik5StatusNeed Info=>Invalid
      Open/ClosedOpen=>Closed
    Wed 25 Nov 2015 12:19:44 AM UTCrik5StatusNone=>Need Info

    Back to the top


    Powered by Savane 3.1-cleanup1