# Mathematica produces seemingly wrong eigenvalues and eigenvectors

I have a Fredholm kernel defined on [-1, 1], which is symmetric under transformation (x, y) -> (-x, -y):

k[x, y] == k[-x, -y]


for any values of x and y, which means that if some function u[x] is an eigenvector of k, then u[-x] is also an eigenvector with exactly the same eigenvalue.

When solving Fredholm eigenvalue equation numerically in Mathematica, the kernel becomes a matrix with exactly the same property:

k[[i, j]] == k[[n + 1 - i, n + 1 - j]]


where n is the number of points. And so the eigenvalues of that matrix must come in pairs. Except that they mostly don't. A few largest eigenvalues do come it pairs, but then something breaks down and Mathematica does not produce pairs anymore. This can be seen by "naked eye" (see below) and so that means that the errors are substantial.

I need to run this for n = 1000. However, the eigenvalues no longer form pairs after the first about 9 pairs even for n = 100, though they hold better for n = 40 for a given example below.

Here is the code that does that along with the check that k[[i, j]] == k[[n + 1 - i, n + 1 - j]] and that the EV decomposition holds. I am using 100 digits of precision but that does not help.

noOfPoints = 100;
maxEigenValues = 40;
precision = 100;
n = 10^4;
m = 4;
e = 0.02;

delta[x_, y_, e_] := 2 * Exp[-(x - y)^2 / e^2]/(e * Sqrt[Pi] * (Erf[(1 - y) / e] + Erf[(1 + y) / e]));
binomial[nn_, kk_] := Gamma[nn + 1.0]/(Gamma[kk + 1.0] * Gamma[nn - kk + 1.0])
entropy[x_, nn_, mm_] := Log[mm^nn * binomial[nn, nn * (x + 1) / 2]] / nn;
rateMultiplier[x_, nn_, mm_] := entropy[0, nn, mm] / entropy[x, nn, mm];

kFunc[x_, y_] := Module[{xp, yp, ep, mp, retVal, np},
xp = SetPrecision[x, precision];
yp = SetPrecision[y, precision];
ep = SetPrecision[e, precision];
mp = SetPrecision[m, precision];
np = SetPrecision[10^4, precision];
retVal = SetPrecision[rateMultiplier[yp, np, mp] * delta[xp, yp, ep], precision];

Return[retVal];
];

domain = {-1, 1};

step = SetPrecision[(domain[[2]] - domain[[1]]) / noOfPoints, precision];
start = SetPrecision[domain[[1]] + step / 2, precision];
midGrid = SetPrecision[Table[start + step * ii, {ii, 0, noOfPoints - 1}], precision];
weights = SetPrecision[Table[step, {ii, 1, noOfPoints}], precision];

k = Table[
SetPrecision[
kFunc[SetPrecision[midGrid[[i]], precision],
SetPrecision[midGrid[[j]], precision]] *
SetPrecision[weights[[i]], precision], precision], {i,
noOfPoints}, {j, noOfPoints}];

Print["k[x, y] - k[-x, -y]"];
kmm = Table[k[[noOfPoints + 1 - ii, noOfPoints + 1 - jj]], {ii, 1, noOfPoints}, {jj, 1, noOfPoints}];
kDiff = k - kmm;
Print["Norm[kDiff] = ", Norm[kDiff]];
ListLinePlot[kDiff, Frame -> True, GridLines -> Automatic, PlotRange -> All]

{val, vec} = Eigensystem[k];
ListPlot[Re[Take[val, Min[maxEigenValues, noOfPoints]]], PlotRange -> All, Frame -> True, GridLines -> Automatic]

Print["Check norm."];
kCheck = k.Transpose[vec] - Transpose[vec].DiagonalMatrix[val];
norm = Norm[kCheck];
Print["norm = ", norm];
ListLinePlot[kCheck, Frame -> True, GridLines -> Automatic, PlotRange -> All]


So, the questions are:

1. Is it possible to use some combination of Method -> ... to make Eigensystem produce a correct result? The documentation for Method is fairly short and my "experiments" did not produce anything useful.
2. If Eigensystem is completely broken, then is it possible to use some other decomposition (which does not have such issue) and then "reconstruct" EV decomposition with correct eigenvalues?
• I think your analysis is wrong, not the eigenvalues. $k(x, y) = k(-x, -y)$ means that the eigenvectors will be either symmetric or antisymmetric; it does not mean "if some function $u(x)$ is an eigenvector of $k$, then $u(-x)$ is also an eigenvector with exactly the same eigenvalue" as you say. For small eigenvalues, the symmetric & antisymmetric eigenvectors come in pairs (negligible symmetry breaking); for larger eigenvalues this pairing breaks up. Nothing wrong with Mathematica's eigenvalue solver. Commented Feb 28, 2022 at 21:29
• @Roman, Oh, no, it is correct. If Integrate[k[x, y] * u[y], {y, -1, 1}] == c * u[x] for any value of x, then c * u[x] == Integrate[k[x, y] * u[y], {y, -1, 1}] == Integrate[k[-x,- y] * u[y], {y, -1, 1}] (because k[x, y] == k[-x, -y]) and then we can apply transformation (x, y) -> (-x, -y) to obtain Integrate[k[x, y] * u[-y], {y, -1, 1}] == c * u[-x]. Call u[-x] = v[x] and we get Integrate[k[x, y] * v[y], {y, -1, 1}] == c * v[x]. Which is exactly what I said. Commented Feb 28, 2022 at 21:49
• At 3.) If the operator k is neither self-adjoint nor normal (it isn't), then how is Mathematica supposed to orthogonalize the eigenbasis? math.stackexchange.com/a/1092999/447001 Commented Feb 28, 2022 at 22:15
• Just two side-notes: 1) Remove all the SetPrecision gymnastics: instead, define e = 0.02100 and replace the 1.0 with 1 in your definition of binomial. 2) You seem to lose performance by reimplementing Binomial in terms of Gamma functions; use Binomial instead (compare Module[{pairs = ReverseSort /@ RandomInteger[{20, 30}, {1000000, 2}]}, {AbsoluteTiming[binomial @@@ pairs;], AbsoluteTiming[Binomial @@@ pairs;], binomial @@@ pairs == Binomial @@@ pairs} ]). Commented Feb 28, 2022 at 22:25
• On a related note, congratulations for casting your frustration into a concise and formal question. This is a very efficient way of study and learning. Commented Mar 1, 2022 at 9:21

I think Roman is right.

Let's denote the partity operator by $$\Pi$$, i.e., $$(\Pi \,u)(x) := u(-x)$$ and the integral operator by $$K$$.

What you proved in your comment: If $$u$$ is an eigenvector of $$K$$, then $$\Pi \, u$$ is also an eigenvector. What you have not shown: that $$\Pi \, u$$ is linearly independent from $$u$$. And in fact, it frequently is not:

See, $$\Pi$$ and $$K$$ commute, i.e., $$K \, \Pi = \Pi \, K$$. Let's define the two operators $$P_\pm := \frac{1}{2}(I \mp \Pi).$$

With $$\Pi \, \Pi = I$$ it is straight-forward to check that these are projectors, i.e., that $$P_+ P_+ = P_+ \quad \text{and} \quad P_- P_- = P_-.$$ Moreover, they satisfy $$P_\pm K = K \, P_\pm, \quad P_+ P_- = P_- P_+ = 0, \quad \text{and} \quad P_+ + P_- = I.$$ In particular this means:

1. The full vector space is a direct sum of the images $$\mathrm{ima}(P_+)$$ and $$\mathrm{ima}(P_-)$$ of $$P_+$$ and $$P_-$$.

2. $$K$$ maps each of $$\mathrm{ima}(P_+)$$ and $$\mathrm{ima}(P_-)$$ to themselves.

Note that $$\mathrm{ima}(P_+) = \mathrm{ker}(P_-)$$ is the subspace of even functions and that $$\mathrm{ima}(P_-) = \mathrm{ker}(P_+)$$ is the subspace of odd functions. That means that we can find an eigenbasis (IIRC, the spectrum of an integral operator can only have 0 as sole accumulation point, so all nonzero spectral values are eigenvalues) solely consisting of even and odd functions. Btw., this does not mean that there won't be any eigenvectors that are neither even nor odd. But it shows that there are many eigenvectors for which $$u$$ and $$\Pi \,u$$ are linearly dependent.

• That's interesting. I plotted 10 pairs of eigenvectors corresponding to the largest 10 eigenvalues for 200 points and, sure enough, the first few ones came in pairs because they were extremely close to zero around x = 0` but then they start to behave very differently around zero and that explains that the eigenvalues can no longer be the same. Thanks for clarification. Commented Feb 28, 2022 at 23:15