Efficiently find all values of parameter such that any of the eigenvalues of a matrix is equal to 1

Question

I wish to find all values for a parameter such that my matrix has an eigenvalue of 1.

Here is an example 16-by-16 matrix with elements depending on the parameter x :

matrix[x_]:= {{8500651/(106043877*(-(34563219/38133806) - x)), 63407530/(1287051647*(104486064/225042547 - x)), 0, 0, 0, 0, 20277381/(169339442*(75256417/41896542 - x)), -(3958441/(26028795502*(152589326/35191063 - x))),34238516/(296067989*(77503175/78528458 - x)),-(53057896/(1058820821*(178917535/75889517 - x))), 0, 0, 0,0,13892954/(196808509*(320486341/58788069 - x)), 23405904/(215304701*(1889864855/236489256 - x))},{53483560/(1085614027*(-(34563219/38133806)-x)),40700493/(163356004*(104486064/225042547 - x)), 0, 0, 0, 0,30006913/(244954743*(75256417/41896542 - x)), -(24764567/(162722462*(152589326/35191063 - x))),-(53057896/(1058820821*(77503175/78528458 - x))), -(87382735/(446966061*(178917535/75889517 - x))), 0, 0, 0, 0, -(42907817/(453704937*(320486341/58788069 - x))), 84963737/(547901533*(1889864855/236489256 - x))},{0, 0, 27658330/(135847301*(61958873/52713692 - x)), -(15114742/(175933119*(206167491/55493486 - x))),-(20277381/(169339442*(-(43728926/153162047) - x))), -(30006913/(244954743*(82298778/75840643 - x))), 0, 0, 0, 0, 53599490/(604925933*(459547987/149753531 - x)), 190081293/(1662771190*(935766167/166848555 - x)),-(13892954/(196808509*(174409325/51756543 - x))), 80315825/(849255191*(202744750/42768967 - x)), 0, 0},{0, 0, -(15114742/(175933119*(61958873/52713692 - x))), 43015171/(201299071*(206167491/55493486 - x)),1011203/(6649182367*(-(43728926/153162047) - x)), 24764567/(162722462*(82298778/75840643 - x)), 0, 0, 0, 0, 190081293/(1662771190*(459547987/149753531 - x)), -(10949853/(262303619*(935766167/166848555 - x))),-(23405904/(215304701*(174409325/51756543 - x))), -(90003800/(580403143*(202744750/42768967 - x))), 0, 0},{0, 0, -(20277381/(169339442*(61958873/52713692 - x))), 995982/(6549096425*(206167491/55493486 - x)),114271858/(1297281265*(-(43728926/153162047) - x)), 13329588/(305900819*(82298778/75840643 - x)), 0, 0, 0, 0, -(19741315/(200069307*(459547987/149753531 - x))), -(19492415/(220036534*(935766167/166848555 - x))),13448559/(144489725*(174409325/51756543 - x)), -(6618957/(182683390*(202744750/42768967 - x))), 0, 0}, {0, 0, -(30006913/(244954743*(61958873/52713692 - x))), 24764567/(162722462*(206167491/55493486 - x)), 24428776/(560616171*(-(43728926/153162047) - x)), 26793265/(204420213*(82298778/75840643 - x)), 0, 0, 0, 0,10786712/(327438171*(459547987/149753531 - x)), -(13724249/(202691483*(935766167/166848555 - x))),-(6618957/(182683390*(174409325/51756543 - x))), -(26501433/(213208117*(202744750/42768967 - x))), 0, 0},{20277381/(169339442*(-(34563219/38133806) - x)), 30006913/(244954743*(104486064/225042547 - x)), 0, 0, 0, 0, 43672883/(160797056*(75256417/41896542 - x)), -(24029253/(358443637*(152589326/35191063 - x))), 19741315/(200069307*(77503175/78528458 - x)), -(17964325/(545319623*(178917535/75889517 - x))), 0, 0, 0, 0, 17162080/(225948483*(320486341/58788069 - x)), 20344659/(116212984*(1889864855/236489256 - x))},{-(995982/(6549096425*(-(34563219/38133806) - x))), -(32693319/(214820528*(104486064/225042547 - x))), 0, 0, 0, 0, -(24029253/(358443637*(75256417/41896542 - x))), 21386598/(109640185*(152589326/35191063 - x)),39606638/(447092233*(77503175/78528458 - x)), 20597410/(304200221*(178917535/75889517 - x)), 0, 0, 0, 0, 52822907/(301735588*(320486341/58788069 - x)), -(11740639/(157158360*(1889864855/236489256 - x)))},{24905741/(215365428*(-(34563219/38133806) - x)), -(53057896/(1058820821*(104486064/225042547 - x))), 0, 0, 0, 0, 25232718/(255722195*(75256417/41896542 - x)), 19492415/(220036534*(152589326/35191063 - x)),20539265/(77909671*(77503175/78528458 - x)), -(23392503/(995388524*(178917535/75889517 - x))), 0, 0, 0, 0, 53720547/(331042420*(320486341/58788069 - x)), 62568183/(552172939*(1889864855/236489256 - x))},{-(50777275/(1013308858*(-(34563219/38133806) - x))), -(87382735/(446966061*(104486064/225042547 - x))), 0, 0, 0, 0, -(14395811/(436994890*(75256417/41896542 - x))), 13724249/(202691483*(152589326/35191063 - x)), -(12404002/(527810181*(77503175/78528458 - x))), 46976883/(200664262*(178917535/75889517 - x)), 0, 0, 0, 0, 15985825/(353717122*(320486341/58788069 - x)), -(132238698/(972444067*(1889864855/236489256 - x)))}, {0, 0, 24073814/(271698003*(61958873/52713692 - x)), 128859719/(1127224173*(206167491/55493486 - x)), -(19741315/(200069307*(-(43728926/153162047) - x))), 17964325/(545319623*(82298778/75840643 - x)), 0, 0, 0, 0, 37555297/(218036843*(459547987/149753531 - x)), 35519067/(482669702*(935766167/166848555 - x)),-(53720547/(331042420*(174409325/51756543 - x))), -(15985825/(353717122*(202744750/42768967 - x))), 0, 0}, {0, 0, 128859719/(1127224173*(61958873/52713692 - x)), -(10949853/(262303619*(206167491/55493486 - x))), -(19492415/(220036534*(-(43728926/153162047) - x))), -(13724249/(202691483*(82298778/75840643 - x))), 0, 0, 0, 0, 35519067/(482669702*(459547987/149753531 - x)), 35641443/(188202652*(935766167/166848555 - x)),-(62568183/(552172939*(174409325/51756543 - x))), 26987636/(198459051*(202744750/42768967 - x)), 0, 0},{0, 0, -(17026893/(241204097*(61958873/52713692 - x))), -(19535558/(179702415*(206167491/55493486 - x))), 26147867/(280929586*(-(43728926/153162047) - x)), -(6618957/(182683390*(82298778/75840643 - x))), 0, 0, 0, 0, -(57288098/(353026759*(459547987/149753531 - x))), -(62568183/(552172939*(935766167/166848555 - x))), 28970920/(162996539*(174409325/51756543 - x)), 30519243/(2860726402*(202744750/42768967 - x)), 0, 0}, {0, 0, 65012549/(687439178*(61958873/52713692 - x)), -(90003800/(580403143*(206167491/55493486 - x))), -(22215436/(613146627*(-(43728926/153162047) - x))), -(23864717/(191995330*(82298778/75840643 - x))), 0, 0, 0, 0, -(15985825/(353717122*(459547987/149753531 - x))), 132238698/(972444067*(935766167/166848555 - x)), 8766113/(821693084*(174409325/51756543 - x)), 136181885/(757167702*(202744750/42768967 - x)), 0, 0}, {17026893/(241204097*(-(34563219/38133806) - x)), -(42907817/(453704937*(104486064/225042547 - x))), 0, 0, 0, 0, 17162080/(225948483*(75256417/41896542 - x)), 36583783/(208974286*(152589326/35191063 - x)), 54543828/(336115729*(77503175/78528458 - x)), 13983559/(309413136*(178917535/75889517 - x)), 0, 0, 0, 0,52664022/(224639591*(320486341/58788069 - x)), 7895426/(230265753*(1889864855/236489256 - x))}, {23405904/(215304701*(-(34563219/38133806) - x)), 90003800/(580403143*(104486064/225042547 - x)), 0, 0, 0, 0, 36583783/(208974286*(75256417/41896542 - x)), -(38839715/(519902359*(152589326/35191063 - x))), 19428068/(171455409*(77503175/78528458 - x)), -(132238698/(972444067*(178917535/75889517 - x))), 0, 0, 0, 0, 7895426/(230265753*(320486341/58788069 - x)), 35132357/(125483278*(1889864855/236489256 - x))}};

I wish to find values of x within a certain range (say -100<x<100) for which any one of the eigenvalues of matrix is 1. The simplest way I can think of to do this is by recognising that if one of the eigenvalues of this matrix is 1, then the determinant of this matrix minus the identity matrix is 0. I find the corresponding x parameters satisfying this by using Reduce:

findParameter = Sort[
                 N[
                   Reduce[Det[matrix[x] - IdentityMatrix[16]] == 0. && -100 <= x <= 100, x, Reals]
  ]
] 

Which gives the correct results:

Out:=
x == -1.0072067712062946 || x == -0.39522236367591385 || x == 0.1697476411232668 || x == 0.7133650412219289 || x == 0.8290303230563018 || x == 1.079792925789695 || x == 1.5468168440854655 || x == 2.1406855069143496 || x == 2.799506930062505 || x == 3.253812273361672 || x == 3.533348897160239 || x == 4.126064806739337 || x == 4.570247581953315 || x == 5.2528098763247995 || x == 5.451459484991106 || x == 7.7272239344443285

My problem is that for larger matrices (for example 100-by-100, of similar sparsity to the example matrix here) this does not work (or it takes far too long), which I suspect is due to the increase in cost of calculating Determinants of larger matrices.

Is there a quicker/more-efficient way to find values for x??

I'm failing at the very first hurdle using Eigenvalues as I can't get FindRoot to work for this 16-by-16 matrix - even knowing the results I'm looking for! For example,

eigenvals[x_] := Eigenvalues[matrix[x]]
FindRoot[eigenvals[x][[1]] == 1., {x, -1}]

does not seem to work.

Answer 1

You can use the Arnoldi-Lanczos algorithm to efficiently find the eigenvalue that is closest to a target value (here, target=1):

closestEVtotarget[x_?NumericQ, target_?NumericQ] :=
  First@Eigenvalues[matrix[N[x]], 1,
  Method -> {"Arnoldi", "Criteria" -> "Magnitude", "Shift" -> target}]

Then it's a matter of plotting and root-finding:

With[{target = 1},
  Plot[closestEVtotarget[x, target], {x, -10, 10}, GridLines -> {None, {target}}]]

With[{target = 1},
  FindRoot[closestEVtotarget[x, target] == target, {x, -1}]]
(* {x -> -1.00721} *)

You can start the root-finder either at hand-picked points (glanced from the plot) or at regularly spaced points:

With[{target = 1},
  Union[Table[x /. FindRoot[closestEVtotarget[x, target] == target, {x, x0}],
    {x0, -2, 9, 1/100}], SameTest -> (Abs[#1 - #2] < 10^-13 &)]]

{-1.00721, -0.395222, 0.169748, 0.713365, 0.82903, 1.07979, 1.54682, 2.14069, 2.79951, 3.25381, 3.53335, 4.12606, 4.57025, 5.25281, 5.45146, 7.72722}

Alternatively we can use the Graphics`Mesh`FindIntersections function (see 199038, 156975, 10475) to get good starting values geometrically from the plot intersections:

With[{target = 1},
  plot = Plot[{target, closestEVtotarget[x, target]}, {x, -10, 10}];
  intersections = Graphics`Mesh`FindIntersections[plot]]

{{-1.00725, 1.}, {-0.958231, 1.}, {-0.395304, 1.}, {-0.351347, 1.}, {0.16972, 1.}, {0.292244, 1.}, {0.713331, 1.}, {0.756746, 1.}, {0.828908, 1.}, {0.941439, 1.}, {1.07977, 1.}, {1.10962, 1.}, {1.54678, 1.}, {1.65261, 1.}, {2.14052, 1.}, {2.222, 1.}, {2.79948, 1.}, {2.94262, 1.}, {3.25374, 1.}, {3.28433, 1.}, {3.53329, 1.}, {3.59908, 1.}, {4.12603, 1.}, {4.21907, 1.}, {4.57021, 1.}, {4.64757, 1.}, {5.2528, 1.}, {5.31701, 1.}, {5.4514, 1.}, {5.52141, 1.}, {7.72721, 1.}, {7.85725, 1.}}

Not all of these are useful: some come from branch jumps. Also, they are not very precise. We refine them with FindRoot:

refined = 
  Union[x /. FindRoot[closestEVtotarget[x, #[[2]]] == #[[2]], {x, #[[1]]}] & /@
    intersections, SameTest -> (Abs[#1 - #2] < 10^-13 &)]

{-1.00721, -0.395222, 0.169748, 0.713365, 0.82903, 1.07979, 1.54682, 2.14069, 2.79951, 3.25381, 3.53335, 4.12606, 4.57025, 5.25281, 5.45146, 7.72722}

Answer 2

for your timing issues you're going to want to run dependencies on x, so that you are doing purely numerical calculations for those larger matrices. The best explanation as to why is given here:

• Issue with computing eigenvalues using Mathematica

I will, however, self promote my answer, which you may find useful here:

• Computing the Eigensystem of a large number of large matrices

So, to overcome your first hurdle, with large matrices, be sure to evaluate your Eigen-related functions on purely numerical elements. And, to overcome your second, you can use a variant of the pure function compiler I provided in the linked answer of mine, which you can also see here

Export[NotebookDirectory[]<>"PureFunctionMatrix.wdx",ToExpression[StringReplace[ToString[UserDefinedMatrixBuildingFunction[a,b,c],InputForm],{"a"->"#1","b"->"#2","c"->"#3"}]<>"&"]];

You can then easily map (or thread) Eigensystem over an imported and precompiled matrix as such

M=Import[NotebookDirectory[]<>"PureFunctionMatrix.wdx"];
eigSet=ParallelMap[Eigensystem[M[#]]&/@xvars]

where xvars is a list of your x-variables of choice. I use Eigensystem here in order to ensure that the outputs of the matrices are orthonormalized, saving me related issues later on, though you would be welcome to use Eigenvalues or Eigenvectors in this case as well, if you have better math skills than I, or are not as lazy. Through my research I have been gathering that a ParallelMap of Eigensystem can negate the integrated parallelization, but I have not yet done extensive testing (again). I will update this entry after such a thing has occurred.

Finally, with a list of your choice Eigenvalues, with each head corresponding to your xvar, you can search this list (pretty efficiently, especially compared to whatever you calculated previous to get the list!) for the eVals that are equal to 1, using this

Table[Pick[Range[Length[eValxlist[xVar]]],eValxlist[xVar],1],{xVar,Length@eValxlist[xVar]}]

which should give you the indexes of the modes in which the corresponding eigenvalues are 1. You can then use these as indexes in order to access your main collation of eVals and eVecs. I prefer this method over the creation of a table, and then the application of a Block[] function to the table of values which searches for what you desire, which is the case in the below linked question from which all of these search methods are inspired, pulled, or referenced from:

Pick elements of largest absolute value

I hope this answer is of some value to you, @Trock, and others, as well! Please let me know if I can clarify anything further, or if I made any egregious (or otherwise!) errors. I know this is not a complete I/o answer, but it will provide a great set of starting points towards solving your problem, and achieving your goals.