Take the 2-minute tour ×
MathOverflow is a question and answer site for professional mathematicians. It's 100% free, no registration required.

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.

share|improve this question
    
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? –  Anthony Quas Jun 19 '12 at 1:03
    
That's the idea (e.g uniformity of the points generated on the circle in $\mathbb{R}^3$) –  David Jun 19 '12 at 7:47
    
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. –  S. Carnahan Jun 20 '12 at 8:36
add comment

marked as duplicate by Igor Rivin, S. Carnahan Jun 20 '12 at 8:35

This question has been asked before and already has an answer. If those answers do not fully address your question, please ask a new question.

5 Answers

up vote 2 down vote accepted

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

share|improve this answer
add comment

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.

share|improve this answer
    
Of course, you should do it w.r.t. a linear subspace if you want the normalization to preserve the subspace. –  Ori Gurel-Gurevich Jun 19 '12 at 4:20
    
How would you normalize w.r.t a linear subspace ? I think here lies the point I am missing. –  David Jun 19 '12 at 8:07
    
You normalize after you project, and the normalized vector stays parallel to the subspace. –  Douglas Zare Jun 19 '12 at 8:14
    
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 ? –  David Jun 19 '12 at 8:22
    
(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). –  Douglas Zare Jun 19 '12 at 9:20
show 5 more comments

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)

share|improve this answer
    
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. –  David Jun 19 '12 at 7:07
add comment

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)?

share|improve this answer
    
it seems difficult to write an algorithm that computes an orthogonal basis of the simplex (and change basis) for any dimension. –  David Jun 19 '12 at 10:01
    
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? –  Goose Jun 19 '12 at 11:32
    
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 ? –  David Jun 19 '12 at 12:12
add comment

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.

share|improve this answer
add comment

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