# Equation solving with magnitude constraints incredibly slow

Somewhat complicated question and sorry in advance if it is presented poorly. I'm running Mathematica 12.0 (haven't updated to 12.1 yet). I'm trying to solve a set of linear complex equations with an additional constraint on the magnitudes of the solutions. What I have tried is very slow, and I'm hoping someone can tell me what I'm doing wrong.

Here is a simplified version of what I tried originally. The real number $$d$$ and the odd prime $$p$$ were defined beforehand. [Edit: I tried to remove all the complicated formulas from the code, but I've added them back in response to a comment.]

p = 3;
d = (p + Sqrt[p^2 + 4*p])/2;
m = 2;
a[n_] := Exp[2*Pi*I*n^2*m/p];
innerProd[g_, h_] := a[g]*a[h]*Conjugate[a[g + h]];
cCubed = 1/Sqrt[p]*Sum[Conjugate[a[g]], {g, 0, p - 1}];
c = cCubed^(1/3);

b[0] = -1/d;
vars = Array[b, p - 1];

equ = {};
For[g = 1, g < p, g++,
AppendTo[equ, b[g] == Conjugate[c*a[g]]/Sqrt[p]*Sum[innerProd[g, h]*b[h], {h, 0, p - 1}]];
]
For[g = 1, g < p, g++,
AppendTo[equ, Abs[b[g]] == 1/Sqrt[p]];
]

Reduce[equ, vars, Complexes]

In the case $$p=3$$, this is a system of two linear equations with two unknowns plus two more equations that specify the magnitude of the solutions. When I tried to run it, Reduce went overnight, used all my memory (16 GB + 32 GB SWAP drive), and did not finish. I believe Solve and NSolve similarly fail.

I managed to solve this issue by splitting the problem into two pieces, as below.

(* same stuff as before *)

b[0] = -1/d;
vars = Array[b, p - 1];

equ = {};
For[g = 1, g < p, g++,
AppendTo[equ, b[g] == Conjugate[c*a[g]]/Sqrt[p]*Sum[innerProd[g, h]*b[h], {h, 0, p - 1}]];
]

solns = Reduce[equ, vars, Complexes] // Simplify;
equ = {};
AppendTo[equ, solns];

For[g = 1, g < p, g++,
AppendTo[equ, Abs[b[g]] == 1/Sqrt[p]];
]

Reduce[equ, vars, Complexes]

What I am doing now is first solving the 2x2 system. That gives me a vector space of solutions in the form of $$b_2 = b_1 + C$$. I am then creating an equation list with that and the magnitude constraints. This runs very quickly and gives me two unique solutions (which is exactly what I want). Does anyone know why Mathematica struggles with the original four equation set when it can do them separately just fine?

Moving on to the case $$p=5$$ (the next easiest case), the first Reduce (a linear system of four equations with four unknowns) runs quickly. The second Reduce, which imposes the magnitude constraints runs overnight without finishing. I had assumed that fixing the magnitudes would not be a difficult part of this computation, so can anyone tell me if I'm doing something horribly wrong?

Just to get out ahead of some possible answers, the next thing I tried was writing each $$b_g$$ as $$\frac{1}{\sqrt{p}}e^{i\theta_g}$$ and solving for $$\theta_g$$ since the magnitudes are known. This way, the equations are no longer linear. They are linear in exponential variables. Still, Mathematica is able to Reduce reasonably quickly. As expected, the answer includes integer constants multiplied by $$2\pi$$. Asking Mathematica to solve the solution it just gave with the additional constraint that $$0 \leq \theta_g < 2\pi$$ for all $$g$$ again goes overnight without finishing. I am thorougly flummoxed as I had assumed limiting to $$[0,2\pi)$$ would be the easy part of this computation. Does anyone know what I'm doing wrong?

Finally, I tried rewriting the linear system as a matrix and using linear algebra functions, but things like LinearSolve give only one solution, and I am looking for all of them. Mathematica documentation says the proper function to use for that is Solve, which brings me back to the original problems I was having.

• Please provide code that readers can run. "function linear in b[i] vars" is not adequate. Sep 30, 2021 at 23:08
• Sorry! I was trying to just ask about the parts I thought were important, but I have added some back. Sep 30, 2021 at 23:37

I'll show something that works quickly on the first example.

First we use ComplexExpand, then replace e.g. Re[b[1]] with re[b[1]] so that it can be treated as a variable (this might not work if we kept the head as the built-in function Re`).

neweq = ComplexExpand[
Apply[Subtract, equ, {1}] // Expand,
{b[1], b[2]}] /. {Re -> re, Im -> im};

Now separate real and imaginary parts.

reeq = ComplexExpand[Re[neweq]];
imeq = ComplexExpand[Im[neweq]];

We can now solve for these over the reals.

Solve[Join[reeq, imeq] == 0, Reals]

(* Out[161]= {{im[b[1]] -> (1/(
9 + 3 Sqrt[21]))(-4 Sqrt[3] -
9 (1/8 (Sqrt[3] - Sqrt[7]) - 1/8 Sqrt[1/3 (6 + 2 Sqrt[21])]) -
3 Sqrt[21] (1/8 (Sqrt[3] - Sqrt[7]) -
1/8 Sqrt[1/3 (6 + 2 Sqrt[21])]) -
1/2 Sqrt[
3] (3 - Sqrt[21] -
6 Sqrt[3] (1/8 (Sqrt[3] - Sqrt[7]) -
1/8 Sqrt[1/3 (6 + 2 Sqrt[21])])) -
1/2 Sqrt[
7] (3 - Sqrt[21] -
6 Sqrt[3] (1/8 (Sqrt[3] - Sqrt[7]) -
1/8 Sqrt[1/3 (6 + 2 Sqrt[21])])) -
1/4 Sqrt[
3] (-3 + Sqrt[21] +
12 Sqrt[3] (1/8 (Sqrt[3] - Sqrt[7]) -
1/8 Sqrt[1/3 (6 + 2 Sqrt[21])])) -
1/4 Sqrt[
7] (-3 + Sqrt[21] +
12 Sqrt[3] (1/8 (Sqrt[3] - Sqrt[7]) -
1/8 Sqrt[1/3 (6 + 2 Sqrt[21])]))),
im[b[2]] ->
1/8 (Sqrt[3] - Sqrt[7]) - 1/8 Sqrt[1/3 (6 + 2 Sqrt[21])],
re[b[1]] ->
1/12 (-3 + Sqrt[21] +
12 Sqrt[3] (1/8 (Sqrt[3] - Sqrt[7]) -
1/8 Sqrt[1/3 (6 + 2 Sqrt[21])])),
re[b[2]] ->
1/6 (3 - Sqrt[21] -
6 Sqrt[3] (1/8 (Sqrt[3] - Sqrt[7]) -
1/8 Sqrt[1/3 (6 + 2 Sqrt[21])]))}, {im[b[1]] -> (1/(
9 + 3 Sqrt[21]))(-4 Sqrt[3] -
9 (1/8 (Sqrt[3] - Sqrt[7]) + 1/8 Sqrt[1/3 (6 + 2 Sqrt[21])]) -
3 Sqrt[21] (1/8 (Sqrt[3] - Sqrt[7]) +
1/8 Sqrt[1/3 (6 + 2 Sqrt[21])]) -
1/2 Sqrt[
3] (3 - Sqrt[21] -
6 Sqrt[3] (1/8 (Sqrt[3] - Sqrt[7]) +
1/8 Sqrt[1/3 (6 + 2 Sqrt[21])])) -
1/2 Sqrt[
7] (3 - Sqrt[21] -
6 Sqrt[3] (1/8 (Sqrt[3] - Sqrt[7]) +
1/8 Sqrt[1/3 (6 + 2 Sqrt[21])])) -
1/4 Sqrt[
3] (-3 + Sqrt[21] +
12 Sqrt[3] (1/8 (Sqrt[3] - Sqrt[7]) +
1/8 Sqrt[1/3 (6 + 2 Sqrt[21])])) -
1/4 Sqrt[
7] (-3 + Sqrt[21] +
12 Sqrt[3] (1/8 (Sqrt[3] - Sqrt[7]) +
1/8 Sqrt[1/3 (6 + 2 Sqrt[21])]))),
im[b[2]] ->
1/8 (Sqrt[3] - Sqrt[7]) + 1/8 Sqrt[1/3 (6 + 2 Sqrt[21])],
re[b[1]] ->
1/12 (-3 + Sqrt[21] +
12 Sqrt[3] (1/8 (Sqrt[3] - Sqrt[7]) +
1/8 Sqrt[1/3 (6 + 2 Sqrt[21])])),
re[b[2]] ->
1/6 (3 - Sqrt[21] -
6 Sqrt[3] (1/8 (Sqrt[3] - Sqrt[7]) +
1/8 Sqrt[1/3 (6 + 2 Sqrt[21])]))}} *)

Reconstructing b[1] and b[2] from this is now straightforward.

--- edit ---

Yes, the case for p=5 does appear to be much more difficult. I changed some of the setup to effectively remove the explicit Conjugate in the resulting equations.

p = 5;
d = (p + Sqrt[p^2 + 4*p])/2;
m = 2;
a[n_] := Exp[2*Pi*I*n^2*m/p];
conja[n_] := Exp[-2*Pi*I*n^2*m/p];
innerProd[g_, h_] := a[g]*a[h]*conja[g + h];
cCubed = 1/Sqrt[p]*Sum[conja[g], {g, 0, p - 1}];
c = FullSimplify[cCubed]^(1/3);
b[0] = -1/d;
vars = Array[b, p - 1];
equ = Join[
Table[-b[g] +
Conjugate[c]*conja[g]/Sqrt[p]*
Sum[innerProd[g, h]*b[h], {h, 0, p - 1}], {g, 1, p - 1}],
Table[Abs[b[g]]^2 - 1/p, {g, 1, p}]];

We can slightly simplify the preprocessing here.

neweq = equ /. {b[n_] -> re[b[n]] + I*im[b[n]]};
reeq = ComplexExpand[Re[neweq]];
imeq = ComplexExpand[Im[neweq]];

But this next step might still hang.

Solve[Join[reeq, imeq] == 0]

Further restricing to Reals will likely make it slower still...

--- end edit ---

• Thank you! This is indeed much quicker, and the answer agrees with what I have for $p=3$. I tried changing $p$ to $5$, and it's having issues. When I run the line to produce neweq, Mathematica complains Internal precision limit $MaxExtraPrecision = 50. Reached while evaluating Im[2e^{-\frac{4i\pi}{5}} + 2e^{\frac{4i\pi}{5}}]. I believe this answer should just be zero, though I tried increasing the MaxExtraPrecision size to 100,000 just in case. I now get the same message, but with 100k instead of 50. Also, the Solve line now just returns itself when$p=5$. Oct 2, 2021 at 9:24 • After above edit: Indeed, Solve still takes a long time for$p=5$. I even tried Solveing the first Table and adding the solution to the second Table, which helped in my original approach. I guess this problem is more difficult than expected. I also tried NSolve. The first time, it told me there were infinitely many solutions and it would intersect with a particular equation, but then gave no solutions. I have not been able to replicate that warning. Since then, it completes but gives an empty set, a result I have reason to be skeptical of. Oct 4, 2021 at 22:18 • After another simplification I can replicate the problem you saw. I'll look more when I get a chance. Oct 4, 2021 at 23:10 • The "infinite solutions" might be a minor bug. I think the real issue is that there are 13 nontrivial equations in 8 unknowns, so the system is over-determined and very possibly just inconsistent. Oct 4, 2021 at 23:35 • I forgot to mention in my last comment, but I also changed the second Table to go$g=1$to$p-1$rather than$p$since$p = 0$in mod$p\$. I did not think the system was inconsistent, but perhaps it is. I will mark the problem as solved for now and look into it some more (outside of Mathematica). Thank you! Oct 5, 2021 at 1:08