eigs routine in octave

Question

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?

What is the actual output you are getting? Possible duplicate of http://scicomp.stackexchange.com/questions/8073/matlab-eigs-wrong-eigenvalues-for-tridiagonal-matrix — Kirill, Jan 20 '15 at 07:34
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. — sebastian_g, Jan 20 '15 at 08:29
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). — Federico Poloni, Jan 20 '15 at 15:28
Dear Federico, thanks for your remark. I modified the code above to check if eigenvalues are real. By the way, although the matrix is not symmetric the eigenvalues are real. This is due to the fact that the analytic solution of the eigenvalue problem has real positive eigenvalues. Regards — sebastian_g, Jan 20 '15 at 18:02
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.) — Bill Greene, Jan 20 '15 at 19:44
@BillGreene You are right, this is likely a bug. Compare for i=1:100; [V, D] = eigs(Lh, 8, 'sr'); fprintf('%f\n', norm(D, 'inf')); end and for i=1:100; D = eigs(Lh, 8, 'sr'); fprintf('%f\n', norm(D, 'inf')); end, I only get bad answers with the [V,D] form. — Kirill, Jan 20 '15 at 23:02
@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. — sebastian_g, Jan 21 '15 at 08:45

score 2 · Answer 1 · answered Apr 06 '18 at 07:32

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.

eigs routine in octave

1 Answers1