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

• To: mathgroup at smc.vnet.net
• 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\$m1u@smc.vnet.net>
• Sender: owner-wri-mathgroup at wolfram.com

```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][]];
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
them.

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
equation.

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,
Prepend[coeffs,lam]]]
]

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:= gvals1 = genEigvals1[A, B, x] Out= {-0.891892, 0}

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

In:= gvals3 = genEigvals3[A,B,x] Out= {-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`"