Equation solving with magnitude constraints incredibly slow

Question

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.

Answer 1

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 ---