Does ordinary "eig" calculate eigenvalues faster than "polyeig"?

I want to compare the computational complexity (computational order) of obtaining eigenvalues (not eigenvectors) of two methods:

• Ordinary eigenvalue calculation, such as "eig" function.
• Polynomial approach, which corresponds to the "polyeig" function.

The main point is that MATLAB consumes less time with "eig" function than "polyeig", however the matrices in the "polyeig" method are smaller in size. Also, there is a line in explanation of "polyeig" which says: " The polyeig function uses the qz factorization to find intermediate results in the computation of generalized eigenvalues. polyeig uses the intermediate results to determine if the eigenvalues are well-determined". Furthermore "eig" uses LU-decomposition in its function, where "polyeig" uses qz decomposition.

My questions are:

1) I think MATLAB's "tic-toc" function is not a proper tool for comparison between these two functions, due to the fact that they uses different methods and "polyeig" checks its answers. Is this suggestion right?

2) For a fair comparison, should I write my own program which calculates eigenvalues in both methods (polynomial and ordinary)?

[Christine Tobler]

Hi Masoud,

Can you tell us more about what problem you are comparing the two functions for? POLYEIG solves a larger set of problems than EIG, and I would certainly expect EIG to be a much better choice than POLYEIG for those problems it supports (simple and generalized eigenvalue problems).

[Masoud]

In fact I want to know when a characteristic equation is more beneficial and requires less computational complexity. I found out that Polyeig uses eig inside its algorithm, which means polyeig can never beat eig, at least in MATLAB platform! In other word, I should compare them in another software, MATLAB is not a proper tool.

[Christine Tobler]

That's correct, all solvers for polynomial eigenvalue problems I am aware of are in some way generalizations of solvers for the simple or generalized eigenvalue problems. They will typically be slower when applied to the same problem.

The solver used in MATLAB calls into EIG directly, with matrices of size n*p, where n is the size of the input matrices, and p is the polynomial degree of the eigenvalue problem.

Since the computational complexity of EIG is O(n^3), the computational complexity of POLYEIG is O(p^3*n^3).

[Walter Roberson]

Masoud comments

Great point

Answer 1

[Pavl M.]

clc
clear all
close all

%format native-bit
format longg
rand('state',1)

%for ProblemComplexty = 2:1:1000

ProblemComplexity = 4;

%[a1,b1] = eig(ones(4,4))

%[a,b,c] = polyeig([ones(1,4);zeros(1,4);zeros(1,4);zeros(1,4)],[zeros(1,4);ones(1,4);zeros(1,4);zeros(1,4)],[zeros(1,4);zeros(1,4);ones(1,4);zeros(1,4)],[zeros(1,4);zeros(1,4);ones(1,4);zeros(1,4)])


A = randn(ProblemComplexity,ProblemComplexity)

A = 4×4

-0.432564811528221 -1.14647135068146 0.327292361408654 -0.588316543014189 -1.6655843782381 1.190915465643 0.174639142820925 2.1831858181971 0.125332306474831 1.1891642016521 -0.186708577681439 -0.136395883086596 0.287676420358549 -0.0376332765933176 0.725790548293303 0.11393131352081

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

B = randn(ProblemComplexity,ProblemComplexity)

B = 4×4

1.06676821135919 0.29441081639264 -0.691775701702287 -1.44096443190102 0.0592814605236053 -1.3361818579378 0.857996672828263 0.571147623658178 -0.095648405483669 0.714324551818952 1.25400142160253 -0.399885577715363 -0.832349463650022 1.62356206444627 -1.59372957644748 0.689997375464345

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

tstart = cputime;

[V,D,W] = eig(A,B)

V = 4×4

1 -1 -0.564096814047548 -0.577094789193665 0.420551114175413 -0.233324018341873 0.208587294472125 0.569423599904986 0.234364513652587 -0.0245396546049334 -1 0.432158863800365 0.960947796142888 -0.53513947825968 -0.220998889485382 -1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

D = 4×4

3.93938105751949 0 0 0 0 -2.89854741781082 0 0 0 0 -0.409481554894135 0 0 0 0 0.470532974588248

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

W = 4×4

0.245304083290022 -0.954105799698367 0.300585105170963 0.442313649434121 1 -1 -0.0419455839846326 -0.0115586114369713 0.57350462982635 -0.228567950330134 -0.196945103378633 1 0.770511810816906 -0.474830861513806 1 0.328909476081192

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

first_method_took_time_sec = cputime - tstart

first_method_took_time_sec = 
        0.0599999999999454

tstart2 = cputime;

[X,e,s] =  polyeig(A,-B)

X = 4×4

1 -1 -0.564096814047548 -0.577094789193665 0.420551114175413 -0.233324018341873 0.208587294472125 0.569423599904986 0.234364513652587 -0.0245396546049334 -1 0.432158863800365 0.960947796142888 -0.53513947825968 -0.220998889485382 -1

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

e = 4×1

3.93938105751949 -2.89854741781082 -0.409481554894135 0.470532974588248

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

s = 4×1

4.03813327484369 5.81249808987839 1.36550981888017 2.18896807862899

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

second_method_took_time_sec = cputime - tstart2

second_method_took_time_sec = 
        0.0300000000002001

speedup = 100*(first_method_took_time_sec - second_method_took_time_sec) / (first_method_took_time_sec)

speedup = 
           49.999999999621

%end for

%Constructed from needing help code by
    %https://independent.academia.edu/PMazniker
    %[phone]
    %https://diag.net/u/u6r3ondjie0w0l8138bafm095b
    %https://github.com/goodengineer
    %https://orcid.org/0000-0001-8184-8166
    %https://willwork781147312.wordpress.com/portfolio/cp/
    %https://www.youtube.com/channel/UCC__7jMOAHak0MVkUFtmO-w
    %https://nanohub.org/members/130066
    %https://pangian.com/user/hiretoserve/
    %https://substack.com/profile/191772642-paul-m

[Masoud Arash]

Thanks for the detailed response.