# Improving my QZ-Algorithm (Include Shifts)

I Need to to solve an generalized Eigenvalue Problem and to compare two Methods (QR and QZ) concerning their convergence rate and execution time.

I started with the Basic QR-Algorithm, implemented in Matlab, tested the results and in the last step i included double-shifts and Deflation strategies.

Now it is time to improve my QZ-Algorithm: I think ist way more complicated, as i did not find any QZ-Algorithm implementations in the Internet. So i grabed the paper of Moler and Stewart and wrote my own implementation. I've been very glad that i did manage to write an algorithm producing the correct results. To have a fair comparison i would like to enhance this algorithm as to include single/double shifts so i can compare the Performance between the QR- and QZ-Algorithm. I think the implicit shift should be more stable, but i want the fastest Solution (less iterations).

You can find the paper here: http://www.dtic.mil/dtic/tr/fulltext/u2/746896.pdf And an extensional paper here: http://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19730018792.pdf

This is my algorithm so far (works correctly, but without shifting, so it is quite slow):

function [lambda, iteration] = QZAlgo(A, B, precision, maxIterations, useHouseholder)
len = length(A);

% 1. Step: Preparation of Matrices
% Transform A to Hessenberg form and B to upper triangular form
if (useHouseholder ~= 0)
[Q, B] = TriangularFormByHouseholder(B);
else
[Q, B] = TriangularFormByGivens(B);
Q = Q';
end
A = Q*A;
[A, B] = HessenbergFormByGivensWhileMaintainingB(A, B);

% 2. Step: QZ-Iteration
for iteration=1:maxIterations
p = 1;
S = A - p*B;

if (useHouseholder ~= 0)
[Q, R] = TriangularFormByHouseholder(S);
else
[Q, R] = TriangularFormByGivens(S);
Q = Q';
end

B1 = Q*B;
[B2, Z] = RepairTriangularFormByGivensFromRight(B1);

B = Q*B*Z;
A = Q*A*Z;

lambda = diag(A) ./ diag(B);

if (CheckForRequiredPrecision(A, B, precision))
break;
end
end
end


The Triangular-Functions are basically QR-Decompositions. The HessenbergFormByGivensWhileMaintainingB-Function creates an QR-Decomposition of A while repairing changes of B with Givens-Transformations of the right side.

The Algorithm works so far, but is, as already mentioned, quite slow.

Do anybody know how i can implement some shifting-strategies to the algorithm?

ANY help, no matter of what Kind and no matter about which shifting strategy, is GREATLY appreciated. I will spend you a crate of beer of your choice, if we can manage to get this working, if you like ;-)

Thanks in Advance!

EDIT: If you Need to have a look in some of the subfuctions, i can provide them here, too.

EDIT2: I tried to have a look at CLAPACK (http://www.netlib.org/clapack/) to figure it out, but the Code is so confusing to me that i had to give it up...

EDIT3: Some of my thoughts about the explicit shift

• Choose shift p: I've Chosen the element in the right bottom Corner of S=A-B as i thought it must be similiar to the procedure of the QR-Algorithm
• Calculate QR-Decomposition of A-p*B
• Calculate B'=QBZ just the same way i did above
• Calculate A=SZ + pB = Q(A-pB) + p*B' But this Procedure did not work. I am clueless at this Point :-(

## 1 Answer

Why don't your consider the "eigs", a built-in function in MATLAB ? Are your matrices nonsymmetric? If symmetric, maybe you can also choose FEAST solver or JDQR.m http://www.staff.science.uu.nl/~vorst102/JDQR.html for solving your generalized eigenvalue problems.

• I do not consider built-in functions of MATLAB, because i am going to implement that algorithm in C. So, i use MATLAB just for testing the correct functionality of the algorithm before i am starting the C implementation. What kind of algorithm is JDQR? If you consider the following generalized eigenvalue Problem $Ax = \lambda Bx$, thent $B$ is symmetric, but $A$ is not! Sep 10 '15 at 11:14
• OK. Maybe, I suggest you look at ARPACK SOFTWARE: caam.rice.edu/software/ARPACK. Besides, it seems that JDQR method can be available for nonsymmetric generalized eigenvalue problems. You can find some references in website staff.science.uu.nl/~vorst102/JDQR.html Sep 10 '15 at 23:27
• It is planned to implement that algorithm on an embedded Hardware so i am not interested in using an rather big library. So, it would like to take my QZ-Algorithm as a Basis for further improvements. Does anyone have an idea on how to implement shifts in the algorithm above? Sep 11 '15 at 14:11
• It is pricy, but for an industrial production setting (vs. research or a one-time use), Mathworks does have the Embedded Coder Oct 10 '15 at 15:07