2
$\begingroup$

I am using octave and observed a problem with the eigs-routine for non symmetric matrices. Using GNU octave version 3.8.1 the code below gives significant difference of eigenvalues although same commands are executed. The difference of eigenvalues is computed by subracting the eig- and eigs-eigenvalues.

% 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');
    if isreal(diag(LambdaEigs))
      LambdaEig=sort(diag(LambdaEig));
      lambdaEigs=sort(diag(LambdaEigs));
    else
      warning('eigenvalues of eigs are not real')
    end
    errmax=max(abs(lambdaEig-lambdaEigs));
    disp(['eigs flag: ' num2str(fl) ', maximum error of eigenvalues: '...
            num2str(errmax)])
    comp=comp+1;
end

Do you have an explanation for this? GNU octave's manual reports that eigs is based on the ARPACK package. I guess ARPACK is widely used for large eigenvalue problems and tested well. But GNU octave's eigs-routine seems not to be reliable for this problem.

I would like to use the eigs routine for eigenvalue computations of larger general eigenvalue problems $A x=\lambda B x$ with $A\in\mathbb{R}^{2000\times2000}$ in GNU octave. Because of the matrix size Krylow subspace methods are reasonably to me and as far as I know eigs is a modified Arnoldi iteration.

Do you have any recommendations how to continue?

Improve this question
$\endgroup$
7
  • 1
    $\begingroup$ What is the actual output you are getting? Possible duplicate of scicomp.stackexchange.com/questions/8073/… $\endgroup$
    – Kirill
    Jan 20, 2015 at 7:34
  • $\begingroup$ Here is the output fl = 0 errmax = 1.5815e+04 fl = 0 errmax = 1.3520e-12 fl = 0 errmax = 2.2737e-12 fl = 0 errmax = 2.7285e-12 fl = 0 errmax = 1.8190e-12 fl = 0 errmax = 1.7906e-12 fl = 0 errmax = 1.5818e+04 fl = 0 errmax = 1.5805e+04 fl = 0 errmax = 2.2737e-12 fl = 0 errmax = 1.3074e-12 fl is the flag of eigs indicating convergence.I also saw the reference you gave. There it says that specifying the number of Lanczos vectors in opts is a solution. But it did not work out. $\endgroup$
    – sebastian_g
    Jan 20, 2015 at 8:29
  • 1
    $\begingroup$ I haven't run the code, but one issue could be that you are calling sort on the eigenvalues, which might be complex in general since your boundary conditions create a nonsymmetric matrix. For complex numbers, sort defaults to ordering by nondecreasing modulus: there could be mismatches if there are eigenvalues with the same modulus (e.g., complex conjugate pairs). $\endgroup$
    – Federico Poloni
    Jan 20, 2015 at 15:28
  • 1
    $\begingroup$ I experimented a bit with your code in octave 3.8.1. The fact that the results vary so widely from trial to trial is very disturbing. I modified your code slightly to print all 8 eigenvalues from both eig and eigs. Sometimes only one or two differ; sometimes all 8 are very different. I would encourage you to submit a bug report to the Octave team. (The MATLAB version of eigs doesn't exhibit this problem.) $\endgroup$
    – Bill Greene
    Jan 20, 2015 at 19:44
  • 2
    $\begingroup$ @BillGreene I already submitted a bug report to the developers. @Kirill It is interesting that bad answers occur only if eigs is called with [V,D]=eigs(...). I will add this to my bug report. Thanks a lot. $\endgroup$
    – sebastian_g
    Jan 21, 2015 at 8:45

1 Answer 1

Reset to default
2
$\begingroup$

As per the submitted bug report, there was a bug in ARPACK < 3.1.5 regarding dneupd, the subroutine that "returns the converged approximations to eigenvalues". In ARPACK 3.2.0, this particular issue has been resolved. Current distributions of octave are usually linked with a much higher version of ARPACK; thus, such behaviour is not expected to happen anymore. But checking that your octave is linked against a new (technically, I want to say stable version of ARPACK) is worth a shot.

I've tested the code above in octave 4.2.2 with ARPACK 3.5.0 installed using MacPorts. The results seem to be stable with a max difference in the eigenvalues between the two around 1E-12, which is, I think, more than satisfactory.

Improve this answer
$\endgroup$

Not the answer you're looking for? Browse other questions tagged or ask your own question.