eigen.pas

delphi矩阵运算、回归分析等数学计算

{ **********************************************************************
  *                           Unit EIGEN.PAS                           *
  *                            Version 2.1d                            *
  *                     (c) J. Debord, March 2004                      *
  **********************************************************************
          Procedures for computing eigenvalues and eigenvectors
  **********************************************************************
  References:
  1) 'Mathematiques et Statistiques' by H. Haut (PSI ed.) : Jacobi
  2) EISPACK (http://www.netlib.org/eispack) : EigenVals, EigenVect
  ********************************************************************** }

unit eigen;

interface

uses
  fmath, matrices;

function Jacobi(A : TMatrix; Lbound, Ubound, MaxIter : Integer;
                Tol : Float; V : TMatrix; Lambda : TVector) : Integer;
{ ----------------------------------------------------------------------
  Eigenvalues and eigenvectors of a symmetric matrix by the iterative
  method of Jacobi
  ----------------------------------------------------------------------
  Input parameters  : A       = matrix
                      Lbound  = index of first matrix element
                      Ubound  = index of last matrix element
                      MaxIter = maximum number of iterations
                      Tol     = required precision
  ----------------------------------------------------------------------
  Output parameters : V      = matrix of eigenvectors (stored by columns)
                      Lambda = eigenvalues in decreasing order
  ----------------------------------------------------------------------
  Possible results  : MAT_OK
                      MAT_NON_CONV
  ----------------------------------------------------------------------
  NB : 1. The eigenvectors are normalized, with their first component > 0
       2. This procedure destroys the original matrix A
  ---------------------------------------------------------------------- }

function EigenVals(A : TMatrix; Lbound, Ubound : Integer;
                   Lambda_Re, Lambda_Im : TVector) : Integer;
{ ----------------------------------------------------------------------
  Eigenvalues of a general square matrix
  ----------------------------------------------------------------------
  Input parameters  : A      = matrix
                      Lbound = index of first matrix element
                      Ubound = index of last matrix element
  ----------------------------------------------------------------------
  Output parameters : Lambda_Re = real part of eigenvalues
                      Lambda_Im = imaginary part of eigenvalues

                      The eigenvalues are unordered, except that complex
                      conjugate pairs appear consecutively with the
                      value having the positive imaginary part first.
  ----------------------------------------------------------------------
  Possible results  :  0 : No error
                      -i : Non-convergence has been encountered during
                           the search for the i-th eigenvalue. The
                           eigenvalues should be correct for indices
                           (i+1)..Ubound.
  ----------------------------------------------------------------------
  NB : This procedure destroys the original matrix A
  ---------------------------------------------------------------------- }

function EigenVect(A : TMatrix; Lbound, Ubound : Integer;
                   Lambda_Re, Lambda_Im : TVector; V : TMatrix) : Integer;
{ ----------------------------------------------------------------------
  Eigenvalues and eigenvectors of a general square matrix
  ----------------------------------------------------------------------
  Input parameters  : A      = matrix
                      Lbound = index of first matrix element
                      Ubound = index of last matrix element
  ----------------------------------------------------------------------
  Output parameters : Lambda_Re = real part of eigenvalues
                      Lambda_Im = imaginary part of eigenvalues
                      V         = matrix of eigenvectors

  If the i-th eigenvalue is real, the i-th column of V contains its
  eigenvector. If the i-th eigenvalue is complex with positive imaginary
  part, the i-th and (i+1)-th columns of V contain the real and imaginary
  parts of its eigenvector. The eigenvectors are unnormalized.
  ----------------------------------------------------------------------
  Possible results  :  0 : No error
                      -i : Non-convergence has been encountered during
                           the search for the i-th eigenvalue. None of
                           the eigenvectors has been found. The
                           eigenvalues should be correct for indices
                           (i+1)..Ubound.
  ----------------------------------------------------------------------
  NB : This procedure destroys the original matrix A
  ---------------------------------------------------------------------- }

procedure DivLargest(V : TVector; Lbound, Ubound : Integer;
                     var Largest : Float);
{ ----------------------------------------------------------------------
  Normalizes an eigenvector V by dividing by the element with the
  largest absolute value
  ---------------------------------------------------------------------- }

function RootPol(Coef : TVector; Deg : Integer;
                 X_Re, X_Im : TVector) : Integer;
{ ----------------------------------------------------------------------
  Real and complex roots of a real polynomial by the method of the
  companion matrix
  ----------------------------------------------------------------------
  Input parameters  : Coef = coefficients of polynomial
                      Deg  = degree of polynomial
  ----------------------------------------------------------------------
  Output parameters : X_Re = real parts of root
                      X_Im = imaginary parts of root
  ----------------------------------------------------------------------
  If no error occurred, the function returns the number of real roots,
  and the roots are sorted in increasing order of their real parts.

  If an error occurred during the search for the i-th root, the function
  returns (-i). The roots should be correct for indices (i+1)..Deg. The
  roots are unordered.
  ---------------------------------------------------------------------- }

implementation

  function Jacobi(A : TMatrix; Lbound, Ubound, MaxIter : Integer;
                  Tol : Float; V : TMatrix; Lambda : TVector) : Integer;
  var
    I, J, K, Im, Jm, Iter : Integer;
    B, C, C2, Na, Nd, P, Q, S, S2, R, T : Float;
  begin
    Iter := 0;
    Na := 0.0;
    Nd := 0.0;
    R := 0.0;

    for I := Lbound to Ubound do
      begin
        V[I,I] := 1.0;
        Nd := Nd + Sqr(A[I,I]);
        if I <> Ubound then
          for J := Succ(I) to Ubound do
            begin
              R := R + Sqr(A[I,J]);
              V[I,J] := 0.0;
              V[J,I] := 0.0;
            end;
      end;

    Na := Nd + 2.0 * R;

    repeat
      R := 0.0;
      for I := Lbound to Pred(Ubound) do
        for J := Succ(I) to Ubound do
          begin
            T := Abs(A[I,J]);
            if T > R then
              begin
                R := T;
                Im := I;
                Jm := J;
              end;
          end;

      B := A[Im,Im] - A[Jm,Jm];

      if B = 0 then
        begin
          C := SQRT2DIV2;
          S := C * Sgn(A[Im,Jm]);
        end
      else
        begin
          P := 2.0 * A[Im,Jm] * Sgn(B);
          Q := Abs(B);
          R := Pythag(P, Q);
          C := Sqrt(0.5 * (1.0 + Q / R));
          S := 0.5 * P / (R * C);
        end;

      for K := Lbound to Ubound do
        begin
          R := V[K,Im];
          V[K,Im] := C * R + S * V[K,Jm];
          V[K,Jm] := C * V[K,Jm] - S * R;
        end;

      if Im <> Lbound then
        for K := Lbound to Pred(Im) do
          begin
            R := A[K,Im];
            A[K,Im] := C * R + S * A[K,Jm];
            A[K,Jm] := C * A[K,Jm] - S * R;
          end;

      if Jm <> Succ(Im) then
        for K := Succ(Im) to Pred(Jm) do
          begin
            R := A[Im,K];
            A[Im,K] := C * R + S * A[K,Jm];
            A[K,Jm] := C * A[K,Jm] - S * R;
          end;

      if Jm <> Ubound then
        for K := Succ(Jm) to Ubound do
          begin
            R := A[Im,K];
            A[Im,K] := C * R + S * A[Jm,K];
            A[Jm,K] := C * A[Jm,K] - S * R;
          end;

      Nd := Nd + 2.0 * Sqr(A[Im,Jm]);

      C2 := Sqr(C);
      S2 := Sqr(S);
      P := 2.0 * S * C * A[Im,Jm];
      R := A[Im,Im];
      A[Im,Im] := C2 * R + S2 * A[Jm,Jm] + P;
      A[Jm,Jm] := S2 * R + C2 * A[Jm,Jm] - P;
      A[Im,Jm] := 0.0;

      Inc(Iter);
      if Iter > MaxIter then
        begin
          Jacobi := MAT_NON_CONV;
          Exit;
        end;
    until Abs(1.0 - Na / Nd) < Tol;

    { The diagonal terms of the transformed matrix are the eigenvalues }
    for I := Lbound to Ubound do
      Lambda[I] := A[I,I];

    { Sort eigenvalues and eigenvectors }
    for I := Lbound to Pred(Ubound) do
      begin
        K := I;
        R := Lambda[I];
        for J := Succ(I) to Ubound do
          if Lambda[J] > R then
            begin
              K := J;
              R := Lambda[J];
            end;
        Swap(Lambda[I], Lambda[K]);
        SwapCols(I, K, V, Lbound, Ubound);
      end;

    { Make sure that the first component of each eigenvector is > 0 }
    for J := Lbound to Ubound do
      if V[Lbound, J] < 0.0 then
        for I := Lbound to Ubound do
          V[I,J] := - V[I,J];

    Jacobi := MAT_OK;
  end;

  procedure Balance(A : TMatrix; Lbound, Ubound : Integer;
                    var I_low, I_igh : Integer; Scale : TVector);
{ ----------------------------------------------------------------------
  This procedure is a translation of the EISPACK procedure Balanc.

  This procedure balances a real matrix and isolates eigenvalues
  whenever possible.

  On input:

    A contains the input matrix to be balanced.

    Lbound, Ubound are the lowest and highest indices
    of the elements of A.

  On output:

    A contains the balanced matrix.

    I_low and I_igh are two integers such that A[i,j]
    is equal to zero if
      (1) i is greater than j and
      (2) j=Lbound,...,I_low-1 or i=I_igh+1,...,Ubound.

    Scale contains information determining the permutations
    and scaling factors used.

    Suppose that the principal submatrix in rows I_low through I_igh
    has been balanced, that P[j] denotes the index interchanged
    with j during the permutation step, and that the elements
    of the diagonal matrix used are denoted by D[i,j].  then
        Scale[j] = P[j],    for j = Lbound,...,I_low-1
                 = D[j,j],      j = I_low,...,I_igh
                 = P[j]         j = I_igh+1,...,Ubound.
    the order in which the interchanges are made is
    Ubound to I_igh+1, then Lbound to I_low-1.

    Note that Lbound is returned for I_igh if I_igh is < Lbound formally
  ---------------------------------------------------------------------- }
  const
    RADIX = 2;  { Base used in floating number representation }

  var
    I, J, M : Integer;
    C, F, G, R, S, B2 : Float;
    Flag, Found, Conv : Boolean;

    procedure Exchange;
    { Row and column exchange }
    var
      I : Integer;
      F : Float;
    begin
      Scale[M] := J;
      if J = M then Exit;

      for I := Lbound to I_igh do
        Swap(A[I,J], A[I,M]);

      for I := I_low to Ubound do
        Swap(A[J,I], A[M,I]);
    end;

  begin
    B2 := RADIX * RADIX;
    I_low := Lbound;
    I_igh := Ubound;

    { Search for rows isolating an eigenvalue and push them down }
    repeat
      J := I_igh;
      repeat
        I := Lbound;
        repeat
          Flag := (I <> J) and (A[J,I] <> 0.0);
          I := I + 1;
        until Flag or (I > I_igh);
        Found := not Flag;
        if Found then
          begin
            M := I_igh;
            Exchange;
            I_igh := I_igh - 1;
          end;
        J := J - 1;
      until Found or (J < Lbound);
    until (not Found) or (I_igh < Lbound);

    if I_igh < Lbound then I_igh := Lbound;
    if I_igh = Lbound then Exit;

    { Search for columns isolating an eigenvalue and push them left }
    repeat
      J := I_low;
      repeat
        I := I_low;
        repeat
          Flag := (I <> J) and (A[I,J] <> 0.0);
          I := I + 1;
        until Flag or (I > I_igh);
        Found := not Flag;
        if Found then
          begin
            M := I_low;
            Exchange;
            I_low := I_low + 1;
          end;
        J := J + 1;
      until Found or (J > I_igh);
    until (not Found);

    { Now balance the submatrix in rows I_low to I_igh }
    for I := I_low to I_igh do
      Scale[I] := 1.0;

    { Iterative loop for norm reduction }
    repeat
      Conv := True;

      for I := I_low to I_igh do
        begin
          C := 0.0;
          R := 0.0;

          for J := I_low to I_igh do
            if J <> I then
              begin
                C := C + Abs(A[J,I]);
                R := R + Abs(A[I,J]);
              end;

          { Guard against zero C or R due to underflow }
          if (C <> 0.0) and (R <> 0.0) then
            begin
              G := R / RADIX;
              F := 1.0;
              S := C + R;

              while C < G do
                begin
                  F := F * RADIX;
                  C := C * B2;
                end;

              G := R * RADIX;

              while C >= G do
                begin
                  F := F / RADIX;
                  C := C / B2;
                end;

              { Now balance }
              if (C + R) / F < 0.95 * S then
                begin
                  G := 1.0 / F;
                  Scale[I] := Scale[I] * F;
                  Conv := False;
                  for J := I_low to Ubound do A[I,J] := A[I,J] * G;
                  for J := Lbound to I_igh do A[J,I] := A[J,I] * F;
                end;
            end;
        end;
    until Conv;
  end;

  procedure ElmHes(A : TMatrix; Lbound, Ubound, I_low, I_igh : Integer;
                   I_int : TIntVector);
{ ----------------------------------------------------------------------
  This procedure is a translation of the EISPACK subroutine Elmhes

  Given a real general matrix, this procedure reduces a submatrix
  situated in rows and columns I_low through I_igh to upper Hessenberg
  form by stabilized elementary similarity transformations.

  On input:

    A contains the input matrix.

    Lbound, Ubound are the lowest and highest indices
    of the elements of A.

    I_low and I_igh are integers determined by the balancing procedure
    Balance. If Balance has not been used, set I_low=Lbound, I_igh=Ubound.

  On output:

    A contains the Hessenberg matrix. The multipliers which were used
    in the reduction are stored in the remaining triangle under the
    Hessenberg matrix.

    I_int contains information on the rows and columns interchanged in
    the reduction. Only elements I_low through I_igh are used.
  ---------------------------------------------------------------------- }
  var
    I, J, M, La, Kp1, Mm1, Mp1 : Integer;
    X, Y : Float;
  begin
    La := I_igh - 1;
    Kp1 := I_low + 1;
    if La < Kp1 then Exit;

    for M := Kp1 to La do
      begin
        Mm1 := M - 1;
        X := 0.0;
        I := M;

        for J := M to I_igh do
          if Abs(A[J,Mm1]) > Abs(X) then
            begin
              X := A[J,Mm1];
              I := J;
            end;

        I_int[M] := I;

        { Interchange rows and columns of A }
        if I <> M then
          begin
            for J := Mm1 to Ubound do
              Swap(A[I,J], A[M,J]);

            for J := Lbound to I_igh do
              Swap(A[J,I], A[J,M]);
          end;

        if X <> 0.0 then
          begin
            Mp1 := M + 1;
            for I := Mp1 to I_igh do
              begin
                Y := A[I,Mm1];
                if Y <> 0.0 then
                  begin
                    Y := Y / X;
                    A[I,Mm1] := Y;
                    for J := M to Ubound do
                      A[I,J] := A[I,J] - Y * A[M,J];