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.

Answer 1

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

Answer 2

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