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][[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
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[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