5
$\begingroup$

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

    enter image description here

  • 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

    enter image description here

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

$\endgroup$
5
  • $\begingroup$ Have you considered FindInstance? $\endgroup$
    – mikado
    Apr 15, 2017 at 18:58
  • 1
    $\begingroup$ The results of your last approach, using Table, certainly are correct. What is your concern with them? $\endgroup$
    – bbgodfrey
    Apr 15, 2017 at 19:00
  • 2
    $\begingroup$ Would RandomPoint do what you're after? $\endgroup$
    – Greg Hurst
    Apr 15, 2017 at 19:53
  • $\begingroup$ @ChipHurst Yes! $\endgroup$
    – Garmekain
    Apr 15, 2017 at 21:29
  • 2
    $\begingroup$ 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. $\endgroup$
    – Igor Rivin
    Apr 15, 2017 at 23:43

2 Answers 2

3
$\begingroup$

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}}
$\endgroup$
1
  • $\begingroup$ It is also biased, two peaks seem to have more points than the other two; check it out $\endgroup$
    – Garmekain
    Apr 15, 2017 at 21:22
1
$\begingroup$

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.