# Optimization for NSolve

I am trying to solve the system of n equations numerically. Basically, I am solving saddle point equations for eigenvalues of hermitean matrices and trying to reproduce Wigner's semicircle law:

n = 10;
eqSystem[n_] :=
Table[x[i] ==
2./n*Sum[If[i != j, (x[i] - x[j])^(-1), 0], {j, 1, n}], {i, 1, n}]

testSystem = eqSystem[n];
vars = Table[x[i], {i, 1, n}];
cond = Table[x[i] < x[j], {i, 1, n}, {j, i + 1, n}] // Flatten;
testSolutions = NSolve[Join[testSystem , cond], vars, Reals];


However, even for n=10 it takes insurmountable amount of time to evaluate the solutions. Is there a way to optimize this code so that it could be evaluated in a reasonable amount of time for large (~1000) n ?

n = 10;
eqSystem[n_] :=
Table[x[i] ==
2./n*Sum[If[i != j, (x[i] - x[j])^(-1), 0], {j, 1, n}], {i, 1, n}]

testSystem = eqSystem[n];
vars = Table[x[i], {i, 1, n}];
cond = Table[x[i] < x[j], {i, 1, n}, {j, i + 1, n}] // Flatten;


You could try turning it into a minimization problem by replacing the heads of each element in testSystem with Subtract (I.e. subtract each lhs from its rhs), and taking the squared sum total of the list of equations. If all the equations are satisfied the squared sum total of lhs-rhs squareResid should be 0:

diffs = Subtract @@@ testSystem;
squareResid = Total[(#[[1]] + #[[2]])^2 & /@ diffs];



On my laptop, this takes about 5 seconds for n=10 and the squared residual is ~6*10^-12 so all equations are pretty close to satisfied. The danger is that these arg values testSolutions don't perfectly satisfy the equations, and depending on the magnitude of the gradient of the objective function the arg values could be far from the true root value, so it's important to verify these results:

({nMinResid, testSolutions} =
NMinimize[Join[{squareResid}, cond], vars]) // AbsoluteTiming
{5.20252, {6.13811*10^-12, {x[1] -> -1.5367, x[2] -> -1.13267,
x[3] -> -0.785613, x[4] -> -0.463586, x[5] -> -0.15335,
x[6] -> 0.15335, x[7] -> 0.463586, x[8] -> 0.785613,
x[9] -> 1.13267, x[10] -> 1.5367}}}


The thing I'd be worried about is the fact that eqSystem has things like (x[i] - x[j])^(-1) in it, so the objective function could explode when x[i]s are close to x[j]s

for n=100, it takes about 100 s and the squared residual is ~6*10^-10 (supressing the output of the arg values):

({nMinResid, testSolutions} =
NMinimize[Join[{squareResid}, cond], vars]); // AbsoluteTiming
(*{104.63, Null}*)
nMinResid
(*5.72012*10^-10*)


I don't have time to run the n=1000 case as of writing (I can update if I do have time later), but this does seem faster than trying to use NSolve

• Thanks for the reply! Did you modify something to evaluate it for n=100? Because it was taking more than half an hour on my laptop to evaluate before I aborted it, and when I tried to evaluate for n=20 a numerous amount of errors occured (I guess the most important one is "Power::infy: Infinite expression 1/0. encountered."). For n=10 I obtain the same result as you did. Commented Apr 16 at 16:00
• I didn't modify anything for n=100. I just ran it again from a fresh kernel. It took 177 s instead of 104 s this time, but the squared error is on the order of 10^-18 (notebook). My \$Version is 14.0.0 for Mac OS X ARM (64-bit) (December 13, 2023) I also ran for n=20 and got no messages and it took 6 s
– ydd
Commented Apr 17 at 0:28
• The Power::infy makes me suspect an x[i]=x[j] somewhere...could it be possible the constraints were changed to x[i] <= x[j] instead of x[i] < x[j]?
– ydd
Commented Apr 17 at 0:32
• No, the contraints are still x[i] < x[j]. Maybe it has something to do with the version of Mathematica, because I was running it on 12.0. I'll get an update and try to evaluate it on the latest version Commented Apr 17 at 11:36
• upd: It works for the 14.0 version Commented Apr 18 at 11:56