MathGroup Archive 1999

[Date Index] [Thread Index] [Author Index]

Search the Archive

Re: Solving A x = Lam B x with A & B singular

  • To: mathgroup at
  • Subject: [mg15474] Re: Solving A x = Lam B x with A & B singular
  • From: Daniel Lichtblau <danl>
  • Date: Mon, 18 Jan 1999 04:22:14 -0500
  • Organization: Wolfram Research, Inc.
  • References: <77f0qi$>
  • Sender: owner-wri-mathgroup at

Below are a few crude ways that may or may not work for relatively
simple examples. I am very doubtful that anything less that full-blown
generalized eigenvalues codes, a la Lapack, will be effective for large
and/or numerically dicey problems.

First method: set up an appropriate interpolating polynomial and find
the roots.

genEigvals1[A_, B_, var_] := Module[{n=Length[A], vals, gendet, rank},
    rank = Length[SingularValues[B][[2]]];
    vals = Map[Det[A-#*B]&, Table[j, {j, 1, rank+1}]];
    gendet = InterpolatingPolynomial[vals, var];
    Chop[Apply[List,Map[Last,NRoots[gendet==0, var]]]]

Second method: Perturb nearby nonsingular examples and then average

genEigvals2[A_, B_, var_] := Module[
        {n=Length[A], dd, eigvals1, eigvals2, eigvals3, eigvals4},
    dd = 10^(-3)*IdentityMatrix[n];
        eigvals1 = Eigenvalues[Inverse[B+dd].A];
        eigvals2 = Eigenvalues[Inverse[B-dd].A];
        eigvals3 = Eigenvalues[Inverse[B+I*dd].A];
        eigvals4 = Eigenvalues[Inverse[B-I*dd].A];
        Chop[(eigvals1 + eigvals2 + eigvals3 + eigvals4)/4]

Third method: Set up and explicitly solve the generalized eigen

genEigvals3[A_, B_, var_] := Module[
        {n=Length[A], j, k,
          coeffs, c, len, vec, lam, tmpvec},
        coeffs = Array[c, n-1];
        For [j=1, j<=n, j++,
                vec[j] = Table [If [k==j, 1, 0], {k,n}]];
        tmpvec = vec[n] + Sum[c[j]*vec[j], {j,n-1}];
        Chop[lam /. Solve[A . tmpvec == lam*B.tmpvec,

For a simple example, I'll use a pair of singular 3x3 matrices.

A = N[{{1,2,3}, {4,5,6}, {7,8,9}}];
B = {{-1, 6, 7}, {8, 2, -5.}, {7, 8, 2}};

In[15]:= gvals1 = genEigvals1[A, B, x] Out[15]= {-0.891892, 0}

In[16]:= gvals2 = genEigvals2[A,B,x] Out[16]= {2.68044, -0.891892, 0}

In[17]:= gvals3 = genEigvals3[A,B,x] Out[17]= {-0.891892, 0}

Note that the second method will give junk even for this example because
moving to a nonsingular example forces a third eigenvalue. Also note
that it is assuming the standard sorting will keep corresponding
eigenvalues in the same position, which might fail if there is a pair
of generalized eigenvalues with nearly equal real parts (if I correctly
remember our ordering of eigenvalues).

The third method will fail if there are eigenvectors with last component
equal to zero. I am fairly certain there are other ways in which these
may give wrong or incomplete results.

Daniel Lichtblau
Wolfram Research

  • Prev by Date: Re: typing starting values in FindMinimum
  • Next by Date: Re: Data fitting with Mathematica 3.0
  • Previous by thread: Solving A x = Lam B x with A & B singular
  • Next by thread: bug in the "Calculus`Limit`"