# Get random solution for equation with $n$ parameters and infinite solutions?

I'm trying to get a random solution for an equation with arbitrary number of parameters which has infinite solutions. Note that the answer must be in n dimensions where n is the number of parameters of the function.

My tries:

• The function can be something like x^2 + y^2 + z^2 == 1, in this case, the solution is a sphere and every point that satisfies that equality is a valid solution. This creates a region that consists of valid solutions, so if we generate a random point in that region we have a valid answer for this problem. That can be done using ybeltukov's answer's RegionDistribution:

implregion = ImplicitRegion[x^2 + y^2 + z^2 + t^2 == 1, {x, y, z, t}];
region = DiscretizeRegion[implregion];
pts = RandomVariate[RegionDistribution[region], 100]


but this raises an error saying DiscretizeRegion only takes arguments with less than 4 dimensions.

• I tried using NSolve to get all the points:

NSolve[x^2 + y^2 + z^2 == 1, {x,y,z}] (* I don't know what I was supposed to see here *)


but this gives a warning message that says:

 Warning: NSolve::infsolns: "Infinite solution set has dimension at least 1.
Returning intersection of solutions with (171802*x-113492*y-121484*z)/178835 == 1."


then it outputs some complex valued solutions, which are not valid (in this context):

{{x->0.252383 - 0.412144I, y->-1.17075 - 0.0795353I, z->-0.0214382 - 0.508549I},
{x->0.252383 + 0.412144I, y->-1.17075 + 0.0795353I, z->-0.0214382 + 0.508549I}}


As you can see, it gets those wanted infinite solutions, and it even selects a random solution, but it is always the same, and contains complex numbers when what I want is one random answer of that solution set which consists only of real values.

• I also tried with Reduce doing:

Reduce[x^2 + y^2 + z^2 == 1]


The result is:

(x == -Sqrt[1 - y^2 - z^2] || x == Sqrt[1 - y^2 - z^2])


and generating some points doing something like

Table[
{RandomChoice[{-1,1}] Sqrt[1 - y^2 - z^2], y, z}/.
{y-> RandomReal[{-1,1}],
z-> RandomReal[{-1,1}]},
{1000}]


gives me incorrect answers because x, y and z are not calculated "at the same time".

Edit:

• Using FindInstance distributes non uniformly the points, I used

FindInstance[x^2 + y^2 + z^2 == 1, {x, y, z}, Reals, 1000]


and the result was

• And the last method distributes non uniformly too, doing

Select[Table[{RandomChoice[{-1, 1}] Sqrt[1 - y^2 - z^2], y, z} /.
{y -> RandomReal[{-1, 1}],
z -> RandomReal[{-1, 1}]},
{1000}],
Im@#[[1]] == 0 &]


gives me

this is because the values are calculated depending on the $yz$ plane and in some places x changes faster than in another places.

• Have you considered FindInstance? – mikado Apr 15 '17 at 18:58
• The results of your last approach, using Table, certainly are correct. What is your concern with them? – bbgodfrey Apr 15 '17 at 19:00
• Would RandomPoint do what you're after? – Chip Hurst Apr 15 '17 at 19:53
• @ChipHurst Yes! – Garmekain Apr 15 '17 at 21:29
• Generating a uniform distribution even on the sphere is highly nontrivial, an ellipsoid is very hard, and you can forget about arbitrary hypersurfaces, so your question is much too general. – Igor Rivin Apr 15 '17 at 23:43

FindInstance works quickly but seems to biased toward points that lie in the xy-plane.

SeedRandom[42];
With[{n = 10},
(FindInstance[x^2 + y^2 + z^2 == 1, {x, y, z}, Reals, n] // N)[[All, All, 2]]]

{{0.790419, -0.108434, 0.602893}, {0.367265, 0.930116, 0.},
{-0.538922, -0.842356, 0.}, {-0.538922, 0.842356, 0.},
{0.55489, -0.831924, 0.}, {0.968064, 0.158416, -0.194311},
{-0.407186, -0.913345, 0.}, {0.978044, 0.0658436, 0.197724},
{0.946108, 0.323852, 0.}, {0.55489, 0.831924, 0.}}


So I recommend

SeedRandom[42];
With[{n = 10},
Module[{y, z},
z = RandomReal[{-1, 1}, n];
y = RandomReal[{-1 + #, 1 - #}] & /@ Abs[z];
MapThread[{Sqrt[1 - #1^2 - #2^2], #1, #2} &, {y, z}]]]

{{0.899547, -0.41092, -0.148189}, {0.972747, 0.0791143, -0.217954},
{0.903025, 0.301654, -0.305861}, {0.882229, 0.46164, -0.0925187},
{0.760152, 0.640032, 0.111927}, {0.701265, 0.574829, -0.421661},
{0.868464, 0.284055, -0.406304}, {0.714105, -0.38114, -0.587185},
{0.922918, 0.161119, -0.349661}, {0.321531, -0.0217459, 0.946649}}

• It is also biased, two peaks seem to have more points than the other two; check it out – Garmekain Apr 15 '17 at 21:22

Region functionality in M11.2+ can handle Ellipsoid objects. For example:

pts = RandomPoint[
RegionBoundary @ Ellipsoid[{0,0,0,0},{1,2,3,4}],
10
]


{{0.702264, -0.506915, -1.00603, -2.29828}, {0.672747, 1.1246, 0.948138, -1.44966}, {-0.543024, 0.140708, -2.00356, -2.01652}, {0.652806, -1.37052, -0.93464, -0.339369}, {-0.38219, -0.873541, -2.1748, 1.48397}, {0.220089, -1.55027, 1.56574, -1.11951}, {-0.38265, 0.746885, 2.49541, -0.596283}, {0.52443, 1.36327, -1.48924, 0.471912}, {-0.212721, 1.43821, 1.50939, 1.71813}, {0.649815, -0.909925, 0.557857, -2.31921}}

Check that they all lie on the Ellipsoid:

Norm/@(pts.DiagonalMatrix[1/Range[4]])


{1., 1., 1., 1., 1., 1., 1., 1., 1., 1.}