How more effectively solve problems with "FindInstance" when number of variables is significantly greater then number of equations

Can you help me with solving the following problem:

There are 16 variables and 3 equations (constraints) on them; is it possible to effectively solve such problems in Mathematica? What is more, there should be a constraint on Abs[(x0 + I y0)]^2+Abs[(v0 + I w0)]^2==1 (for other variables similarly).

b0 = {{1}, {0}};
b1 = {{0}, {1}};

these are 4 qubits

psi = (x0 + I y0) b0 + (v0 + I w0) b1;
psi = (x1 + I y1) b0 + (v1 + I w1) b1;
psi = (x2 + I y2) b0 + (v2 + I w2) b1;
psi = (x3 + I y3) b0 + (v3 + I w3) b1;

Example steps for psi

1. we define a density matrix in the following way:

ro = Refine[psi.ConjugateTranspose[psi], {x0, y0, v0, w0} > 0]

2. this is an element of the POVM, 2 times 2 matrix

R = 1/2ro;

3. sum of these matrices should give an identity matrix

sum = Sum[R[i], {i, 0, 3}]

The above requirement allows us to write the following equations(?)

Do you know is it possible to solve such equations with NSolve or find some solutions with FindInstance?

FindInstance[-2 + x0^2 + x1^2 + x2^2 + x3^2 + y0^2 + y1^2 + y2^2 +
y3^2 == 0 && v0 x0 + v1 x1 + v2 x2 + v3 x3 + w0 y0 + w1 y1 + w2 y2 + w3 y3 ==
0 && w0 x0 + w1 x1 + w2 x2 + w3 x3 - v0 y0 - v1 y1 - v2 y2 -
v3 y3 == 0, {x0, y0, v0, w0, x1, y1, v1, w1, x2, y2, v2, w2, x3,
y3, v3, w3}, Reals, 1] // N

I also tried to solve this problem with different parameterization of states:

psi = Cos[θ0/2] b0 + Exp[I*ϕ0] Sin[θ0/2] b1;

Similar steps allow for deriving the following set of equations

FindInstance[
Cos[θ0] + Cos[θ1] + Cos[θ2] + Cos[θ3] ==
0 && Cos[ϕ0] Sin[θ0] + Cos[ϕ1] Sin[θ1] +
Cos[ϕ2] Sin[θ2] + Cos[ϕ3] Sin[θ3] == 0 &&
Sin[θ0] Sin[ϕ0] + Sin[θ1] Sin[ϕ1] +
Sin[θ2] Sin[ϕ2] + Sin[θ3] Sin[ϕ3] ==
0, {θ0, θ1, θ2, θ3, ϕ0, ϕ1, ϕ2, ϕ3}, Reals, 1] // N

This approach doesn't allow me to find anything. Maybe you know how one can solve such problems more effectively?

Thanks for any help with this. J.

• Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory tour now, 2) 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, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise. Dec 17 '21 at 14:42
• FindInstance tries to find exact (algebraic) roots. Do you look for exact solutions?
– Acus
Dec 17 '21 at 15:08

If number of variables is greater than the number of equations on might look for a minimum solution (similar to PseudoInverse for linear problems):

Similar to @Akku14's clever answer try

mini = NMinimize[{    vars1 . vars1, Thread[eqs1 == 0]}, vars1 ]
(*{2., {x0 -> 0.476697, y0 -> -0.857636, v0 -> -1.03538*10^-30,
w0 -> -1.82424*10^-30, x1 -> 0.293517, y1 -> -0.446633,
v1 -> 5.42342*10^-31, w1 -> 9.86076*10^-32, x2 -> 0.68577,
y2 -> -0.295276, v2 -> 5.42342*10^-31, w2 -> 6.40949*10^-31,
x3 -> -0.0942965, y3 -> -0.430379, v3 -> 4.19082*10^-31,
w3 -> -2.46519*10^-31}}*)
• Very clever (+1). Then FindInstance can easily work with the remaining equation eqs1 /. Thread[ Evaluate[Cases[mini[] // Chop, Rule[aa_, 0] -> aa]] -> 0] Dec 17 '21 at 19:52

You can use NMinimize.

vars1 = {x0, y0, v0, w0, x1, y1, v1, w1, x2, y2, v2, w2, x3, y3, v3,
w3};
(eqs1 = Subtract @@@
List @@ (-2 + x0^2 + x1^2 + x2^2 + x3^2 + y0^2 + y1^2 + y2^2 +
y3^2 == 0 &&
v0 x0 + v1 x1 + v2 x2 + v3 x3 + w0 y0 + w1 y1 + w2 y2 +
w3 y3 == 0 &&
w0 x0 + w1 x1 + w2 x2 + w3 x3 - v0 y0 - v1 y1 - v2 y2 -
v3 y3 == 0));

nmin1 = NMinimize[{eqs1.eqs1, Thread[eqs1 == 0]}, vars1,
WorkingPrecision -> 25]

(*   {0, {x0 -> -0.2688698621225149789707399,
y0 -> -0.6778829684493148842221356,
v0 -> -0.4599201839580140286276144,
w0 -> -0.7244698615595609677369415,
x1 -> 0.7764510656415366048951560,
y1 -> -0.4594336779699087692804516,
v1 -> -0.1147213848446481707078258,
w1 -> -0.1394270088263146193311412,
x2 -> -0.004129737804174234957776942,
y2 -> 0.1133188085181842516758908,
v2 -> 0.4228333046093424543366667,
w2 -> 0.4599659401886634735379034,
x3 -> 0.7759590155226886589669734,
y3 -> -0.1981350944948345426075089,
v3 -> -0.6731961155868888303000809,
w3 -> 0.5942934315955578978673677}}   *)

eqs1 /. nmin1[]

(*   {0.*10^-25, 0.*10^-25, 0.*10^-25}   *)

With different methods or starting vaues for vars, you get other solutions.

• Nice answer. Sometimes one might be interested in a minimum solution NMinimize[{ vars1 . vars1, Thread[eqs1 == 0]}, vars1, WorkingPrecision -> 25] Dec 17 '21 at 17:50

FindInstance appears to begin guesses with integer values of 1 or 0. As a result, the first of the three equations yields complex solutions for the remaining variables. So, replace 2 by 200 in the first equation and then renormalize the result

FindInstance[-200 + x0^2 + x1^2 + x2^2 + x3^2 + y0^2 + y1^2 + y2^2 + y3^2 == 0
&& v0 x0 + v1 x1 + v2 x2 + v3 x3 + w0 y0 + w1 y1 + w2 y2 + w3 y3 == 0
&& w0 x0 + w1 x1 + w2 x2 + w3 x3 - v0 y0 - v1 y1 - v2 y2 - v3 y3 == 0,
{x0, y0, v0, w0, x1, y1, v1, w1, x2, y2, v2, w2, x3,y3, v3, w3}] /.
Rule[a_, b_] -> Rule[a, b/10]

(* {x0 -> -(1/10), y0 -> -(1/10), v0 -> 1/10, w0 -> 1/10, x1 -> 1/10,
y1 -> 0, v1 -> 0, w1 -> 1/10, x2 -> 1/10, y2 -> 1/30 (-2 - 4 Sqrt),
v2 -> 0, w2 -> -(1/10), x3 -> 1/20 (-2 + 1/3 (2 + 4 Sqrt)),
y3 -> 1/20 (2 + 1/3 (-2 - 4 Sqrt)), v3 -> -(1/10), w3 -> 1/10} *)

Test the correctness by

Simplify[List @@ (-2 + x0^2 + x1^2 + x2^2 + x3^2 + y0^2 + y1^2 + y2^2 + y3^2 == 0
&& v0 x0 + v1 x1 + v2 x2 + v3 x3 + w0 y0 + w1 y1 + w2 y2 + w3 y3 == 0
&& w0 x0 + w1 x1 + w2 x2 + w3 x3 - v0 y0 - v1 y1 - v2 y2 - v3 y3 == 0) /. %]

(* {True, True, True} *)

Obtain other solutions by selecting different values of the FindInstance option, RandomSeeding.

The same method by @bbgodfrey.

At first we deal with the second and third equations since they are the homogeneous equation.(Furthormore we can Reduce it )

Here we also use some compact notations.

x = {x0, x1, x2, x3};
y = {y0, y1, y2, y3};
v = {v0, v1, v2, v3};
w = {w0, w1, w2, w3};
vars = Flatten[{x, y, v, w}];
ins = FindInstance[{v . x + w . y == 0, w . x - v . y == 0,
And @@ UnequalTo /@ vars}, vars, Reals,1];
norm = x . x + y . y /. ins[];
{x, y, v, w} = Sqrt[2/norm]*{x, y, v, w} /. ins[];
sol = Thread[vars -> Flatten[{x, y, v, w}]]
-2 + x0^2 + x1^2 + x2^2 + x3^2 + y0^2 + y1^2 + y2^2 + y3^2 == 0 &&
v0 x0 + v1 x1 + v2 x2 + v3 x3 + w0 y0 + w1 y1 + w2 y2 + w3 y3 == 0 &&
w0 x0 + w1 x1 + w2 x2 + w3 x3 - v0 y0 - v1 y1 - v2 y2 - v3 y3 ==
0 /. sol

{x0 -> -Sqrt[(2/11)], x1 -> -Sqrt[(2/11)], x2 -> -2 Sqrt[2/11], x3 -> -Sqrt[(2/11)], y0 -> -Sqrt[(2/11)], y1 -> -Sqrt[(2/11)], y2 -> -Sqrt[(2/11)], y3 -> -Sqrt[(2/11)], v0 -> -Sqrt[(2/11)], v1 -> -Sqrt[(2/11)], v2 -> -Sqrt[(2/11)], v3 -> -Sqrt[(2/11)], w0 -> -Sqrt[(2/11)], w1 -> -Sqrt[(2/11)], w2 -> -9 Sqrt[2/11], w3 -> 16 Sqrt[2/11]}