# VERY slow to solve a nonlinear equation involving matrix eigenvalues

I am trying to solve a problem similar to the following:

f[k_, a_] := (
Matrix = {{0.05, (a*k)^2}, {(a*k)^2, -0.05}};

{eps, psi} = Eigensystem[Matrix];

{eps,psi} = {eps[[#]], psi[[#]]} &@
Ordering[eps];
Return[psi[[2]]];
);

FindRoot[f[k, 3] == 5, {k, 0.1}];

Clear["Global*"];


Basically I want to plot a band structure and find out at which parameter a given band crosses the zero energy.

This simple example works great for either FindRoot or NSolve. However, if I replace "Matrix" with a more complicated matrix, like a 300-by-300 Hermitian matrix, neither FindRoot nor NSolve can give me anything. They both stuck and just froze.

For your information, the function f[k, a] in the real problem takes about 0.2 seconds to evaluate each time, which is spent partly on constructing the matrix and partly on the diagonalization process. Moreover, the problem structure is actually simple, as it should always have a unique real solution.

At the same time, both Python and Matlab can get me an answer pretty quickly. Moreover, even if I implement the root-searching algorithm (Bisection or secant method) in Mathematica myself, I can also obtain a solution pretty quickly.

Do you have an idea why this root-searching process cannot be solved efficiently with the built-in functions in Mathematica?

• (1) FindRoot might handle this better if you clear f and then define as f[k_?NumberQ,a_?NumberQ]:=.... (2) Rather than compute all eigenpairs, if you just want the smallest maybe do Eigensystem[matrix,1] Jan 3 '16 at 21:36
• Maybe I'm thinking too far ahead, but is the matrix coming from differential equations? If so, have you looked at DEigenvalues and NDEigenvalues? Jan 4 '16 at 8:16
• @xzczd Thanks for the comments. No, it is not coming from a differential equation. It is just an ordinary Hamiltonian that describes the bandstructure of a certain system.
– Xiao
Jan 4 '16 at 16:43
• @DanielLichtblau Thanks for the reply! Actually your suggestions solved the problem for me. In particular, I've never thought of using NumberQ to constraint the function argument. Is it true that when carrying out purely numerical works in Mathematica, it is better to add constraints like NumberQ to speed up things?
– Xiao
Jan 4 '16 at 16:57
• Regarding _?NumberQ usage and speedup potential, it varies (alas). When symbolic preprocessing e.g. to compute analytic Jacobians or Hessians is the cause of slowness then it helps considerably. My guess is that was the case with the FindRoot as originally defined, hence the restriction was helpful for speed. Jan 4 '16 at 17:06

Here's a solution, based in part on Daniel Lichtblau's comments. Instead of returning the eigenvector psi[[2]], I thought f should return the eigenvalue.

f[k_?NumericQ, a_?NumericQ] :=
Eigenvalues[
{{0.05, (a*k)^2},
{(a*k)^2, -0.05}},
-1];

FindRoot[f[k, 3] == 5, {k, 0.1}] // AbsoluteTiming

(*  {0.007437, {k -> 0.745337}}  *)


Now comparing with the OP's code, also replacing psi[[2]] by eps[[2]], I see no difference in timing (V11.2):

ClearAll[f];
f[k_, a_] := (
Matrix = {{0.05, (a*k)^2}, {(a*k)^2, -0.05}};
{eps, psi} = Eigensystem[Matrix];
{eps, psi} = {eps[[#]], psi[[#]]} &@Ordering[eps];
Return[eps[[2]]];);

FindRoot[f[k, 3] == 5, {k, 0.1}] // AbsoluteTiming

(*  {0.007368, {k -> 0.745337}}  *)
`

Maybe something has changed in the last couple of years.