5
$\begingroup$

I have positions and radii of $n$ $n$-dimensional hyperspheres and want to find their intersection points efficiently. A very-straight-forward solution seems quite reliable:

Timing@With[{d = 50},
  With[{
    p = RandomReal[{-1, 1}, d],
    s = RandomReal[{-1, 1}, {d, d}]},
   Sort[EuclideanDistance[p, #] & /@
     (Array[x, d] /.
       Quiet@Solve[
         Element[Array[x, d], Sphere[#, EuclideanDistance[p, #]]] & /@ s,
         Array[x, d]])]]]

{69.002, {4.84698*10^-14, 0.858948}}

The problem here is that 50-dimensional solution already takes about 70 seconds on my system, and I'd want to find solutions when $n \gg 100$. FindMinimum or NMinimize of sums of squares of sphere surface distances might work to an extent, but these find only one of the solutions. Are there good alternatives to the straight-forward Solve?

Improve this question
$\endgroup$
5
  • $\begingroup$ Are you looking for all-pairs solutions? Also what exactly should a solution contain? Orthogonal vector defining the hyperplane of intersection plus center and radius of the intersecting hypercircle? $\endgroup$
    – Daniel Lichtblau
    Commented Oct 16, 2022 at 16:27
  • $\begingroup$ @DanielLichtblau I'm honestly wondering if I'm naive, but don't almost all intersections of $n$ $n$-spheres result at most two point solutions? In practice I'm looking for those, not pairwise hypercircles (if I stated my intent incompletely). $\endgroup$
    – kirma
    Commented Oct 16, 2022 at 17:21
  • $\begingroup$ I misunderstood; I thought you were seeking pairwise sphere intersection sets. You are quite correct, and apologies for the noise. $\endgroup$
    – Daniel Lichtblau
    Commented Oct 17, 2022 at 14:04
  • 1
    $\begingroup$ Also: Now that I actually understand the question, the linear algebra approach is the method of choice for this. $\endgroup$
    – Daniel Lichtblau
    Commented Oct 17, 2022 at 14:07
  • $\begingroup$ @DanielLichtblau No problem! I also thought that linear algebra solution should be the way to go, but frankly, as this is not what I do for a living I just felt that somebody here probably can quickly come about with a highly performing alternative to the trivial problem statement in Wolfram Language. :) $\endgroup$
    – kirma
    Commented Oct 17, 2022 at 16:00

2 Answers 2

Reset to default
6
$\begingroup$

If pts is a $d \times d$ matrix whose rows are the $d$ center points, and rs is a list of $d$ positive numbers containing the radii, then this will calculate the two intersection points (assuming they intersect at all...):

find[pts_,rs_]:=With[{A=MapThread[hyperplane[First[pts],First[rs],#1,#2]&,{Rest[pts],Rest[rs]}]},
    Module[{s},With[{c=LinearSolve[A[[;;,1]],A[[;;,2]]]+s*First[NullSpace[A[[;;,1]]]]},
        c/.Solve[(c-First[pts]).(c-First[pts])==First[rs]^2,s]]]];
hyperplane[pA_,rA_,pB_,rB_]:={pB-pA,(rA^2-Dot[pA,pA]-rB^2+Dot[pB,pB])/2};

Example

(* generate spheres *)
d=500;
SeedRandom[1];
p0=RandomReal[{-1,1},{d}];
pts=RandomReal[{-1,1},{d,d}];
rs=Map[Norm[#-p0]&,pts];

(* example *)
AbsoluteTiming[Map[Norm[#-p0]&,find[pts,rs]]]
(* {0.148522,{6.9176*10^-12,0.201555}} *)

Comment: Given two spheres $A$ and $B$ with centers $p_A$, $p_B$ and radii $r_A$, $r_B$ then a point $x$ lies on both spheres if $$ \|x-p_A\|^2 = r_A^2 \qquad \|x-p_B\|^2 = r_B^2 $$ Subtracting the two equations gives the linear equation $$ x \cdot (p_B-p_A) = (r_A^2 - \|p_A\|^2 - r_B^2 + \|p_B\|^2)/2 $$ as a necessary condition where the dot means the Dot-product. So this restricts $x$ to a hyperplane. The code above computes the equations for $d-1$ such hyperplanes (intersection of sphere 1 with each of the $d-1$ other spheres), then computes their intersection which generically is a line (see LinearSolve and NullSpace), then computes those two points on the line that lie on sphere 1 (using Solve).

Improve this answer
$\endgroup$
0
2
$\begingroup$

Well, speeding it up significantly (but not as much as @user293787...) was easier than I thought:

Timing@With[{d = 500},
  With[{
    p = RandomReal[{-1, 1}, d],
    s = RandomReal[{-1, 1}, {d, d}]},
   Sort[EuclideanDistance[p, #] & /@
     ({#, 2 RegionNearest[
             Hyperplane[Cross @@ (# - First@s & /@ Rest@s), First@s], #] - #} &@
       (Array[x, d] /.
         Last@FindMinimum[
           Sum[
            ((SquaredEuclideanDistance[Array[x, d], c] /. Abs -> Identity) -
               SquaredEuclideanDistance[p, c])^2, {c, s}],
           Array[x, d], MaxIterations -> 1000]))]]]

(* {54.6429, {7.00921*10^-10, 0.0262783}} *)

So, with FindMinimum and increased MaxIterations one can typically accomplish meaningful results for 500-dimensional case in comparable time to the original 50-dimension case. The original example takes now about 0.1 seconds instead of over a minute.

The second point is found simply by projecting the found solution through the hyperplane spanned by the centres of spheres to the opposite side of the hyperplane.

Improve this answer
$\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.