# How can I speed up my code finding the parameter which yields zero determinant for increasingly larger matrices?

I build a square matrix with matrix elements dependent on the parameter x. I need to find all values of x within a reasonable range (say -10<x<10) such that the determinant of this matrix is 0. My method takes far too long for larger matrices.

First, the matrix elements are built from the numerical eigenvalues and eigenvectors of a previous matrix as follows:

(* In what follows a, b and c are numeric - example values could be a=0.5, b=0.5 and c=-0.2*)

eigens[a_, b_] := eigens[a, b] =
Block[{},
m = 10;
ham = Table[i^2*KroneckerDelta[i, j], {i, -m, m},{j, -m, m}] +
Table[a*KroneckerDelta[i, j + 1], {i, -m, m}, {j, -m, m}] +
Table[a*KroneckerDelta[i, j - 1], {i, -m, m}, {j, -m, m}] +
Table[b*KroneckerDelta[i, j + 2], {i, -m, m}, {j, -m, m}] +
Table[b*KroneckerDelta[i, j - 2], {i, -m, m}, {j, -m, m}];

SortBy[Transpose[Eigensystem[ham, Method -> "Banded"]], First]
] (* Eigenvalues and functions later used in calculating matrix elements *)

val1[a_, b_, n_] := val1[a, b, n] = eigens[a, b][[n + 1]][] (* Eigenvalues *)
func1[a_, b_, n_] := func1[a, b, n] = Chop[(1./Sqrt[N[2*Pi]])*eigens[a, b][[n + 1]][]].Table[Exp[I*j*\[CurlyPhi]], {j, -m, m}]; (* Sum of exponentials where coefficients come from the eigenvector *)
val2[a_, b_, n_] := val2[a, b, n] = val1[-a, -b, n] (* Defined as eigenvalues for negative a and b *)
func2[a_, b_, n_] := func2[a, b, n] = func1[-a, -b, n] (* Same but for the function *)


My matrix elements are then calculated using val1, val2, func1, and func2 via the following:

me[(a_)?NumericQ, (b_)?NumericQ, c_, x_, k_, l_, p_, q_] := Block[{},
prefactor = (-2*Pi*c)/(val1[a, b, p] + val2[a, b, q] - x);
numInt = NIntegrate[
Chop[func1[a, b, p]*func2[a, b, q]*Conjugate[func1[a, b, k]]*Conjugate[func2[a, b, l]]],
{\[CurlyPhi], -Pi, Pi},
AccuracyGoal -> 10,
Method -> {"GlobalAdaptive", Method -> "GaussKronrodRule","SymbolicProcessing" -> 0}
];

prefactor*numInt
]//Chop;


Is there a better way of performing this numerical integral?

Finally I build my square matrix of size (max+1)^2 via:

matrix[a_, b_, c_, x_, max_] := Block[{},
g = Table[me[a, b, c, x, i, f, j, k], {i, 0, max},{f, 0, max}, {j, 0, max}, {k, 0, max}];

h = Rationalize[ArrayReshape[g, {(max + 1)^2, (max + 1)^2}]-IdentityMatrix[(max + 1)^2], 0]
];


and I find the values of x (and sort them in ascending order) such that the determinant is 0 using Reduce:

findParameter[a_, b_, c_, max_] := findParameter[a, b, c, max] =
Sort[N[Reduce[Det[matrix[a, b, c, x, max]] == 0. && -10 <= x <= 10, x, Reals]]]


Is there a more efficient way than using Reduce?

My problem is that for larger and larger matrices this method takes way way too long. For example while findParameter takes 0.4s for max=2, it rapidly scales up and when max=4, findParameter takes 30s . Ideally I would like to evaluate this problem for max=6 quickly, as further calculations currently seems to be taking hours/days!

• Is there a way to go from an answer in n dimensions to an answer in n+1 dimensions? If so, this might allow you to leverage previous answers. – bill s Apr 5 at 15:36
• Are a, b, c meant to stay symbolic all the time? Then there is probably no chance to get it working in reasonable time. Otherwise, please provide typically values for a, b, and c. Guessing parameters is a pain in the neck. – Henrik Schumacher Apr 5 at 15:54
• Moreover, it might help little to use ham = SparseArray[ { Band[{1, 1}] -> Range[-N[m], N[m], 1]^2, Band[{1, 2}] -> N[a], Band[{2, 1}] -> N[a], Band[{1, 3}] -> N[b], Band[{3, 1}] -> N[b] }, {2 m + 1, 2 m + 1}, 0. ] for generating the matrices. – Henrik Schumacher Apr 5 at 15:59
• @HenrikSchumacher a, b, and c are always numeric and something like between 0 and 1 (and c is negative). An example would be a=0.5, b=0.5, and c=-0.2 - I'll edit the question appropriately thanks. – Trock Apr 5 at 17:57