picking a random solution to linear equations

This question is for a specific problem but should be of more general interest. I have two points known in 3D space say [a1 b1 c1] and [a2 b2 c2]. I need to pick two additional (unknown) points say [x1 y1 z1] and [x2 y2 z2] that satisfy the following three equations:

euclidean_distance([a1 b1 c1],[x1 y1 z1]) = 1 and euclidean_distance([a2 b2 c2],[x2 y2 z2]) = 1 and euclidean_distance([x1 y1 z1],[x2 y2 z2]) = 1

so I suppose three equations and six unknowns (the coordinates of the two unknown points). Here is where im at,

if ( euclidean_distance([a2 b2 c2],[a1 b1 c1]) > 3 ) there are 0 solutions. if ( euclidean_distance([a2 b2 c2],[a1 b1 c1]) = 3 ) there is one easy to find solution if ( euclidean_distance([a2 b2 c2],[a1 b1 c1]) < 3 ) there are infinite solution and I want to known one of them at random; preferably evenly picked from the solution space.

I know matlab can do this because matlab can do anything, Thanks

[John D'Errico]

These are not linear equations at all. In fact, they are quadratic, because the points in question lie on the surfaces of spheres, and intersections thereof. As it turns out, the set representing the intersection of two spheres does lie in a plane, but it is still a circular set in that plane, so again, not truly linear.

Answer 1

[Alfonso Nieto-Castanon]

I like your Matlab can do anything remark. The evenly space condition is a bit hard to achieve, but here some thoughts anyway.

I would go about this by parameterizing things with respect to d, the distance between your point X1 (x1 y1 z1) and A2 (a2 b2 c2), and of course D, the distance between your point A1 (a1 b1 c1) and A2 (a2 b2 c2).

D is fixed for your problem. For any value of d, the "number" of points X1 that satisfy your condition (1) and are at a distance d of A2 is proportional to sqrt(4-(D^2-d^2+1)^2/D^2). For each of these points X1, the "number" of points X2 satisfying conditions (2) and (3) is proportional to sqrt(4-d^2). That means that combined, the "number" of solutions for each value of d is proportional to p(d) = sqrt((4-(D^2-d^2+1)^2/D^2)*(4-d^2)).

With that at hand (assuming for a moment that I did not made any mistakes in these equations; please double-check), you could do the following:

1. pick a value of d at random following a distribution with pdf proportional to p(d)
2. select a point X1 at random from the circle intersecting the 1-sphere around A1 and the d-sphere around A2
3. select a point X2 at random from the circle intersecting the 1-sphere around X1 and the 1-sphere around A2.

That, I believe, should give you the desired "random evenly in solution space" solution.

To do the first step (get a random value d with the proper distribution), you could integrate p(d) and do that analytically, but for simplicity I would probably do something like:

 e = abs(D-1):1e-5:min(2,D+1);   % range of possible d values
 p = sqrt(max(0,4-(D^2-e.^2+1).^2/D^2).*max(0,4-e.^2)); % pdf
 P = cumsum(p)/sum(p);           % CDF
 d = e(find(P>rand,1,'first'));  % random d value
 if isempty(d), d=e(1); end      % takes care of D=0 and D=3 extreme cases

The second and third steps (picking X1 and X2 from the value of d) are just a matter of selecting a couple of points randomly from two circles, so I will leave that to you as an exercise. (with Matlab you can do anything as well!)

[Alfonso Nieto-Castanon]

Solution space example: Black dots are A1 and A2, and surfaces represent the color-coded density of X1 and X2 solutions. This example is for D=2.5 (assuming that all my random thoughts above were correct, but I would not put much weight on that; I just like pretty pictures and my simulations are taking a long time today...)

Answer 2

[Igor]

Thanks! very thorough. In the meantime you got me thinking the the right direction and I came up with a simple (albeit wasteful) solution

A = [0 0 1]; B = [0 0 0];

num_test_points = 300; %for larger values the solution is more and more accurate

S1 = random_points_on_sphere(num_test_points,A); %unit test points around A

S2 = random_points_on_sphere(num_test_points,B); %likewise

distances = abs(abs(pdist2(S1,S2))-1); %should be relatively even in the solution space [r,c]=find(distances==min(min(distances)));

P1 = S1(r,:); %point one

P2 = S2(c,:); %point two

solution = [A;P1;P2;B];

plot3(solution(:,1),solution(:,2),solution(:,3)) hold on

[Alfonso Nieto-Castanon]

Yes, that is very clever and perfectly valid.

The only somewhat inconvenient issue there is that you are not guaranteed to find an approximate solution within a certain tolerance of a true solution, rather you get a bound on the likelihood of this happening (and you control that bound with num_test_points). Depending on your application you might want to repeat your procedure until min(min(distances)) is below a predefined threshold, and that should give you a more meaningful way to control the accuracy/tolerance of the resulting solutions.

Another alternative would be to, instead of generating your S1 and S2 points randomly, use a predefined resolution-level triangulation of the sphere surface (e.g. icosahedral partitioning) which would naturally place a hard bound on the tolerance of your resulting approximate solutions.

Just some random thoughts, very interesting question by the way! (just curious, where does this problem arise / come from?)

[Igor]

Im re-doing monte carlo simulations of a polymer from a recnet paper and this is one of the monte carlo moves i.e. replacement of part of the self avoiding walk with a new segment of three bonds. I thought about just doing axle rotation to "scramble" three of the existing bonds but im trying to be as true to the algorithm as possible.