# 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 ;-)

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 :-(

• 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! – Roland Sep 10 '15 at 11:14