2
$\begingroup$

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 :-(
$\endgroup$
-1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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! $\endgroup$ – Roland Sep 10 '15 at 11:14
  • $\begingroup$ 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 $\endgroup$ – Hsien-Ming Ku Sep 10 '15 at 23:27
  • $\begingroup$ 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? $\endgroup$ – Roland Sep 11 '15 at 14:11
  • $\begingroup$ It is pricy, but for an industrial production setting (vs. research or a one-time use), Mathworks does have the Embedded Coder $\endgroup$ – GeoMatt22 Oct 10 '15 at 15:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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