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

Question

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.

The results of your last approach, using Table, certainly are correct. What is your concern with them? — bbgodfrey, Apr 15 '17 at 19:00
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

m_goldberg · Answer 1 · 2017-04-15T20:01:20.073

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

score 1 · Answer 2 · answered Sep 28 '17 at 00:45

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

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

2 Answers2