0
$\begingroup$

Possible Duplicate:
Efficiently sampling points uniformly from the surface of an n-sphere

I'm trying to generate random points of a (n-2)-sphere on a n-hyperplane so basically the intersection of a (n-1)-sphere and a n-hyperplane.

The hyperplane here is the plane associated with a (n-1)-simplex of equation: $x_1 + x_2 + ... + x_n = 1$

I know the radius $R$ of the n-sphere and I know the center $a$ which is in the simplex. $a = (a_1, a_2, ..., a_n)$ and $a_1+a_2+...+a_3=1$

So the points I want to generate are solution of :

$x_1 + x_2 + ... + x_n = 1$ and $(x_1-a_1)^2 + (x_2-a_2)^2 + ... + (x_n-a_n)^2 = R^2$

I absolutely need to generate points uniformly and in Cartesian coordinates.

I have been looking into it for a few days now and cannot seem to find a good way to do it.

For the moment I did it on $\mathbb{R}^3$ which is equivalent to generate points of a circle on a plane. To do so, I projected the space on the 2-simplex (3D->2D) and I generated points of a circle using :

$(u,v) = \frac{R}{\sqrt{X_1^2 + X_2^2}}(X_1,X_2)$ where $X_1, X_2$ are independent gaussians.

(this algorithm allows you to generate points on a n-sphere too)

Finally, I used an "inverse-projection" of the simplex to the space (after some calculus...) :

$z = \sqrt{2/3} \cdot v $

$y = \frac{u}{\sqrt{2}} - 1/2 \cdot z$

$x=1-y-z$

Here is a graph to give you an exemple of what I just presented in 3 dimensions. The standard 2-simplex is the area represented in cyan.

http://www.freeimagehosting.net/t/aksf6.jpg (sorry cannot post image tags...)

Of course, those formulas are easily derived in 3 dimensions but not as easily in $\mathbb{R}^n$.

I feel like I am looking at the problem the wrong way so any idea is welcome.

Thanks.

Improve this question
$\endgroup$
3
  • $\begingroup$ when you say uniformly, you mean uniformly with respect to the $(n-2)$-dimensional Lebesgue measure living on the intersection of the hyperplane and the sphere? $\endgroup$
    – Anthony Quas
    Commented Jun 19, 2012 at 1:03
  • $\begingroup$ That's the idea (e.g uniformity of the points generated on the circle in $\mathbb{R}^3$) $\endgroup$
    – David
    Commented Jun 19, 2012 at 7:47
  • $\begingroup$ The problem of rotating and translating points on a sphere is outside the scope of this forum. You may want to ask at one of the other websites listed in the FAQ. $\endgroup$
    – S. Carnahan ♦
    Commented Jun 20, 2012 at 8:36

5 Answers 5

Reset to default
3
$\begingroup$

I'll answer using matrix vector notation, which generalizes the problem to multiple linear constraints. We wish to generate a point $\mathbf{x}$ which is in the intersection of a hypersphere and an affine space, which we can write like this:

$$\|\mathbf{x}-\mathbf{x}_0\|=r$$ $$\mathbf{Ax}=\mathbf{b}$$

Then we can generate $\mathbf{x}$ with this procedure:

$$\mathbf{z} \sim \mathcal{N}(\mathbf{0},\mathbf{I}) $$ $$\mathbf{z} := \mathbf{z} - \mathbf{A}^\dagger\mathbf{Az}$$ $$\mathbf{y} := \mathbf{A}^\dagger (\mathbf{b}-\mathbf{Ax}_0)$$ $$\mathbf{x} := \mathbf{x}_0 + \mathbf{y} + \frac{\mathbf{z}}{\|\mathbf{z}\|}(r - \|\mathbf{y}\|)$$

The idea here is that the point $\mathbf{x}_0 + \mathbf{y}$ is the center of the intersection, and $\mathbf{z}/\|\mathbf{z}\|$ is uniformly distributed on the unit sphere of the null space of $\mathbf{A}$.

Note, we must have $r \ge \|\mathbf{y}\|$, otherwise the intersection is empty. Also, we assume $\mathbf{A}$ is full rank so that the pseudoinverse satisfies $\mathbf{A}^\dagger = \mathbf{A}^T(\mathbf{A A}^T)^{-1}$. In your case $\mathbf{A} = \mathbf{1}^T$ so $\mathbf{A}^\dagger = n^{-1}\mathbf{1}$.

Improve this answer
$\endgroup$
5
$\begingroup$

You can always generate points uniformly on an n-sphere by sampling n+1 independent gaussians and then normalizing the resulting vector. So you ought to be able to do this by choosing n independent gaussians, projecting orthogonally to the plane, and then normalizing the vector.

Improve this answer
$\endgroup$
10
  • $\begingroup$ Of course, you should do it w.r.t. a linear subspace if you want the normalization to preserve the subspace. $\endgroup$
    – Ori Gurel-Gurevich
    Commented Jun 19, 2012 at 4:20
  • $\begingroup$ How would you normalize w.r.t a linear subspace ? I think here lies the point I am missing. $\endgroup$
    – David
    Commented Jun 19, 2012 at 8:07
  • $\begingroup$ You normalize after you project, and the normalized vector stays parallel to the subspace. $\endgroup$
    – Douglas Zare
    Commented Jun 19, 2012 at 8:14
  • $\begingroup$ I understood the principle but how to apply it (in my case) ? + another question comes to my mind, how do you manage the radius then ? $\endgroup$
    – David
    Commented Jun 19, 2012 at 8:22
  • $\begingroup$ (1) Generate a rotationally symmetric Gaussian distribution. (2) Project to the hyperplane through the origin parallel to your hyperplane. (3) Normalize to get a uniform distribution on a sphere in the hyperplane. (4) Translate back to your hyperplane. Which step is the problem? And if you need a sphere of a different radius you rescale in step (3). $\endgroup$
    – Douglas Zare
    Commented Jun 19, 2012 at 9:20
0
$\begingroup$

See Efficiently sampling points uniformly from the surface of an n-sphere (I have voted to close this question as an exact duplicate of the above, btw)

Improve this answer
$\endgroup$
1
  • $\begingroup$ My problem is not about generating points on a n-sphere. The algorithm in the question you posted is the same as the one I exposed (that's what I meant by "this algorithm allows you to generate points on a n-sphere too"). The problem is more about "rotating" and "translating" a generated sphere to fit on/in the hyperplane. $\endgroup$
    – David
    Commented Jun 19, 2012 at 7:07
0
$\begingroup$

Can't you just find an orthogonal basis for your simplex, generate random points for a (n-2) sphere in cartesian coordinates, and then change basis (and translate)?

Improve this answer
$\endgroup$
3
  • $\begingroup$ it seems difficult to write an algorithm that computes an orthogonal basis of the simplex (and change basis) for any dimension. $\endgroup$
    – David
    Commented Jun 19, 2012 at 10:01
  • $\begingroup$ Why? The hyperplane is parallel to $x_1+x_2+...+x_n = 0 $. So if your sphere is on that plane, a translation will put it on the plane that you want, no? $\endgroup$
    – Goose
    Commented Jun 19, 2012 at 11:32
  • $\begingroup$ Then you need to rotate your generated sphere to fit on the plane $x_1+x_2+...+x_n=0$, and how would you do that ? $\endgroup$
    – David
    Commented Jun 19, 2012 at 12:12
0
$\begingroup$

Ok sorry, I have been not doing geometry for too long !

I found the solution I wanted with the following algorithm which is I think not an exact replica of what Jason and Douglas were proposing but it seems to work :

(for 1 point - just sample m observations if you need more points)

1) Sample 1 observation of n independent gaussians $(a_1,a_2,...,a_n)$ You get a point $X= (x_1,x_2,...,x_n)$

2) Project it orthogonally to the hyperplane (in my case $x_1+x_2+...+x_n=1$) to get the projection $Y=(y_1,y_2,...,y_n)$

3) Normalize $Y$ w.r.t. the projection of the origin.

In my case the projection of the origin is $(1/n, 1/n, ..., 1/n)$

So to normalize I used : $ \displaystyle \frac{Y}{\sqrt{(y_1-1/n)^2 + (y_2-1/n)^2+...+(y_n-1/n)^2}}$

4) Project the normalization one last time orthognally to your hyperplane to get a point on the surface of the unit (n-2)-sphere centered in $(1/n,...,1/n)$

If you want a radius $R$ and a specific center $C$ on the hyperplane you just need to scale with $R$ and add $C$ in 3).

Any comments on that method ? In 3D, when I normalize in 3) I don't get a parallel circle as mention above by Douglas but a sort of cylinder which projected give me the circle. Any explanation ?

Thanks a lot for your help anyway.

Improve this answer
$\endgroup$

Not the answer you're looking for? Browse other questions tagged .