# Dealing with 6 VERY large equations. Am I doing this correctly, and can I expect a solution?

I'm new to Mathematica, but so far, everything seems straight forward and intuitive, which is really nice. However, I've come to use this software since I have to show a rather nasty 3x3 matrix is unitary. Initially, I set up a system to solve ConjugateTranspose[M]*M == IdentityMatrix[3], but that was taking too long, so I decided to verify that the columns are orthonormal, or rather for which values of my parameters solve such a system. Here is the first cell of the code I am dealing with.

    (* Defines the values of the final transfer matrix and the row
vectors corresponding to such a system *)

Rp11 = (U11*U22*beta + U21*U12*T21*V12)/(alpha*beta);
Rp12 = U12/beta;
Rp13 = (T11*T22*V21*U12)/(alpha*beta);
Rp21 = U21/beta;
Rp22 = alpha/beta;
Rp23 = (T12*beta + T11*T22*V12)/beta;
Rp31 = (V11*V22*T21*U21)/(alpha*beta);
Rp32 = (V21*beta + V11*V22*T21)/beta;
Rp33 = (T11*T22*V11*V22)/(alpha*beta);

q1 = {Rp11, Rp12, Rp13};
q2 = {Rp21, Rp22, Rp23};
q3 = {Rp31, Rp32, Rp33};


So, not too bad. The Rpij are the elements of the matrix which I have to show is unitary. The real problem comes in when I take the dot products of the q vectors and try to make them orthonormal, since I have to define the Tij, Uij, and Vij in the following way:

(* Defines the variables in their full glory. *)

alpha = T11*U22*V11;
beta = 1 - T21*V12;
B1 = eta1\[Conjugate]*tau1\[Conjugate] - Exp[I*theta1];
B2 = eta2\[Conjugate]*tau2\[Conjugate] - Exp[I*theta2];
B3 = eta3\[Conjugate]*tau3\[Conjugate] - Exp[I*theta3];
theta1 = phi11 + phi21;
theta2 = phi12 + phi22;
theta3 = phi13 + phi32;
T11 = (eta1\[Conjugate] - tau1*Exp[I*theta1])/B1;
T12 = (gamma1\[Conjugate]*kappa1*Exp[I*phi11])/B1;
T21 = (gamma1*kappa1\[Conjugate]*Exp[I*phi21])/B1;
T22 = (tau1\[Conjugate] - eta1*Exp[I*theta1])/B1;
U11 = (eta2\[Conjugate] - tau2*Exp[I*theta2])/B2;
U12 = (gamma2\[Conjugate]*kappa2*Exp[I*phi12])/B2;
U21 = (gamma2*kappa2\[Conjugate]*Exp[I*phi22])/B2;
U22 = (tau2\[Conjugate] - eta2*Exp[I*theta2])/B2;
V11 = (eta3\[Conjugate] - tau3*Exp[I*theta3])/B3;
V12 = (gamma3\[Conjugate]*kappa3*Exp[I*phi13])/B3;
V21 = (gamma3*kappa3\[Conjugate]*Exp[I*phi23])/B3;
V22 = (tau3\[Conjugate] - eta3*Exp[I*theta3])/B3;
eta1 = er1 + I*ei1;
tau1 = tr1 + I*ti1;
gamma1 = gr1 + I*gi1;
kappa1 = kr1 + I*ki1;
eta2 = er2 + I*ei2;
tau2 = tr2 + I*ti2;
gamma2 = gr2 + I*gi2;
kappa2 = kr2 + I*ki2;
eta3 = er3 + I*ei3;
tau3 = tr3 + I*ti3;
gamma3 = gr3 + I*gi3;
kappa3 = kr3 + I*ki3;


The etas, taus, gammas, and kappas also follow the relations:

    Abs[eta1]^2 + Abs[gamma1]^2 == 1;
Abs[eta2]^2 + Abs[gamma2]^2 == 1;
Abs[eta3]^2 + Abs[gamma3]^2 == 1;
Abs[tau1]^2 + Abs[kappa1]^2 == 1;
Abs[tau2]^2 + Abs[kappa2]^2 == 1;
Abs[tau3]^2 + Abs[kappa3]^2 == 1;


So altogether, I've counted 15 equations and 18 unknowns, so the system should be underspecified, and so sets of values for 3 of the parameters should work. So, at first I tried using Solve, but I read somewhere on this site that Solve expects there to be an exact solution, but since the system is under-specified, I assume there will be sets of values which will allow the inner products between the q vectors to be orthonormal. Then I tried using Reduce and specifying that I want the numbers defining the gammas, taus, kappas, and etas to be real using the Element command, like so:

    Q11 = q1.q1\[Conjugate];
Q12 = q1.q2\[Conjugate];
Q13 = q1.q3\[Conjugate];
Q22 = q2.q2\[Conjugate];
Q23 = q2.q3\[Conjugate];
Q33 = q3.q3\[Conjugate];

Reduce[Q11==1 && Q12==0 && Q13==0 && Q22==1 && Q23==0 && Q33==1 && Element[{er1, ei1, tr1, ti1, gr1, gi1, kr1, ki1, er2, ei2, tr2, ti2, gr2, gi2,kr2, ki2, er3, ei3, tr3, ti3, gr3, gi3, kr3, ki3}, Reals],{er1,tr1,kr1},Reals]


I think I let this run for about 2 days before I gave up on it. I've also tried just simplifying the elements of the Rpij matrix. I'm not sure if I am doing this in a way which 'makes sense' to mathematica. Have I set this up correctly? And if so, should I expect to find a solution?

• Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Read the faq! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! Jul 24, 2015 at 18:08
• @Anon Just some thoughts. I copied your equations and noticed an asymmetry in one spot. You use phi32 where (knowing nothing about the equations) using symmetry I would expect the term to be phi23. I changed that and attempted to run the Reduce that you have at the end. I got an error about using Reals for the Domain (I think it wants Complexes). I see 30 variables, 24 variables named er1,ei1,...kr3,ki3 and six variables phi11, ... phi23. Don't you need to add the absolute value of the eta and gamma relationships in the Reduce? Jul 25, 2015 at 3:19
• Hey Jack. That phi32 should be a phi 23, thanks for the spot! Also, that's a good point. I thought I could define those relations in a separate cell, and then only include the two relations in Reduce that I did. Although, from my limited experimentation, it seems Mathematica needs me to be VERY explicit in specifying what it wants me to do, so I probably shouldn't assume it understands that I want to reduce the equations along with the square norm values in a previous cell. Thank you for the response!!
– Anon
Jul 30, 2015 at 1:27

This appears to give you one approximate solution in a few seconds

sol=NMinimize[
Norm[Abs[eta1]^2 + Abs[gamma1]^2 - 1] + Norm[Abs[eta2]^2 + Abs[gamma2]^2 - 1] +
Norm[Abs[eta3]^2 + Abs[gamma3]^2 - 1] + Norm[Abs[tau1]^2 + Abs[kappa1]^2 - 1] +
Norm[Abs[tau2]^2 + Abs[kappa2]^2 - 1] + Norm[Abs[tau3]^2 + Abs[kappa3]^2 -1 ] +
Norm[Dot[q1, Conjugate[q1]] - 1] + Norm[Dot[q1, Conjugate[q2]]] +
Norm[Dot[q1, Conjugate[q3]]] + Norm[Dot[q2, Conjugate[q2]] - 1] +
Norm[Dot[q2, Conjugate[q3]]] + Norm[Dot[q3, Conjugate[q3]] - 1],
{phi11, phi12, phi13, phi21, phi22, phi23, phi32, er1, ei1, er2, ei2,
er3, ei3, gr1, gi1, gr2, gi2, gr3, gi3, kr1, ki1, kr2, ki2, kr3,
ki3, tr1, ti1, tr2, ti2, tr3, ti3}]


{0.00181976, {phi11 -> 0.851471, phi12 -> 0.175396, phi13 -> 0.421107, phi21 -> -0.408143,
phi22 -> -0.645298, phi23 -> -0.304258, phi32 -> -0.304633, er1 -> -0.79449,
ei1 -> -0.539209, er2 -> -0.832498, ei2 -> 0.554027, er3 -> -0.682537, ei3 -> -0.719161,
gr1 -> -0.029367, gi1 -> -0.277806, gr2 -> -0.00118513, gi2 -> -0.000426324,
gr3 -> 0.129221, gi3 -> 0.0159043, kr1 -> -0.0287896, ki1 -> -0.174254, kr2 -> 0.463902,
ki2 -> 0.451902, kr3 -> -0.224306, ki3 -> 0.375873, tr1 -> -0.482919, ti1 -> -0.857669,
tr2 -> 0.118905, ti2 -> 0.752623, tr3 -> 0.301435, ti3 -> -0.847079}}


Then

{q1, q2, q3}.Conjugate[{q1, q2, q3}] /. sol[[2]]


yields

{{1.00179+2.20106*10^-7 I,    -0.000630293+0.000303027 I, 0.0000315854-0.000011917 I},
{0.000629258+0.000303853 I,   0.998984+0.000452428 I,   -0.0442241-0.00871245 I},
{0.000038223+6.23644*10^-6 I, 0.0442018-0.00870806 I,    0.998984-0.000452648 I}}


which gives some support for the unitary conjecture, at least up to the precision of the result.

You can give NMinimize a WorkingPrecision option to make it work with more digits, or try adjusting MaxIterations or AccuracyGoal, but a few quick experiments didn't show a substantial decrease in the minimum.

You can try to check this approximate solution and see if it is good enough for what you need to do.

If I try @Jack's possible replacement of phi32 with phi23 I get a different solution with a somewhat larger sum of normed errors. If you can resolve that question and any other possible issues then perhaps this method can rapidly get you a useful solution.

• Hey Bill, thanks for taking the time to address my issue; I really appreciate it! Directly to the matter, @Jack was correct in his replacement, and I didn't pick up on that permutation until now. I've fixed it. But what exactly did you do here? It looks like you found the smallest such values of the parameters I've listed above to solve this system. Is there any way I'd be able to determine the set of points for which Unitarity will hold? I mean that since we have more variables than equations, in principle we ought to be able to solve for most in terms of a few, giving surfaces as solutions.
– Anon
Jul 30, 2015 at 1:34
• @Anon: You were having problems trying to find a way to solve when you had fewer equations than variables and with your method not finishing in two days. I just showed one method of getting around both of those in three seconds. I found ONE set of parameters that approximately satisfied your equations. But since the result doesn't rapidly converge to almost exactly the identity matrix I worry there may still be more small errors hidden in it. If your actual question is to find the equations or surfaces that solve this then you may have a much bigger task than I can help you with.
– Bill
Jul 30, 2015 at 8:12
• How do you figure that I have more equations than unknowns? In my view, only the theta equations, square norm equations, and inner products carry any real information, since the others are just defined for convenience (makes substitution easier). Dealing with just the magnitudes of the etas, etc. I find that I have 15 equations and 18 unknowns. If I include the parametrization of the etas, etc. into reals and imaginaries, I have 27 equations and 30 unknowns. So I see that I ought to be able to solve for everything in terms of 3 variables. Thanks again for the help!
– Anon
Jul 30, 2015 at 13:33