2
$\begingroup$

Suppose $A$ is an $n$-sphere of unit radius and $B$ is a set of $x\in\mathbb{R}^{n+1}$ satisfying $\langle a, x\rangle=b$, how do I sample from $A\cap B$?

There's R code for $n=2$ given here but I'm interested in the approach for general $n$.

There are a few heuristic methods implemented here, but unclear how if they respect the original density of points being uniform on the sphere.

$\endgroup$
2
  • 3
    $\begingroup$ Hm. Isn't that intersection just an $n-1$-sphere? So is this not "simply" a case of sampling from that lower-dimensional sphere and mapping back to $n+1$-space? $\endgroup$ Oct 17, 2023 at 16:09
  • 2
    $\begingroup$ Agreed with @Stephan. There's nothing new here. A very simple solution is to use Gram-Schmidt or just a generic rotation to make $a$ the unit vector in the $n+1$st coordinate and generate points on a sphere in the first $n$ coordinates. $\endgroup$
    – whuber
    Oct 17, 2023 at 18:50

2 Answers 2

4
$\begingroup$

There's R code for $n=2$ given here but I'm interested in the approach for general $n$.

I wrote that R code in that way in order to avoid the simulation being about the sub-question of the intersection of a plane and a n-sphere, being a (n-1)-sphere, and to be make the code resemble the problem description instead of using tricks that make the code run faster (yet more difficult to relate with the original question).

The little Easter eggs in that code with the plots like plot(xs[2,],xs[3,]) demonstrate how that way of sampling actually ends up in a (n-1)-sphere (which could have been simulated more directly from the start). I would not recommend this way of rejection sampling in general.


To perform the sampling more quickly and more generally in higher dimensions, you can sample a (n-1)-sphere (of size $\sqrt{1-(b/|a|)^2}$), then transform onto a basis perpendicular to the vector $a$.

$\endgroup$
3
$\begingroup$

Idea

I will adjust your notation a bit to distinguish between vectors and scalars. Use $\vec{a}$ as the fixed vector of unit length in $\mathbb{R}^{n+1}$ and $0 \leq b < 1$ as the constant.

The intersection is an $n-1$ dimensional sphere of radius $\sqrt{1 - b^2}$.

So you could:

  • Use Gram-Schmidt to find an orthonormal basis $\vec{u}_1, \vec{u}_2, \dots \vec{u}_n$ of the hyperplane $\langle \vec{a}, \vec{x} \rangle = 0$.
  • Draw random samples $\vec{c}$ from an $n$ dimensional multivariate normal distribution.
  • Divide those by their length, and then multiply by $\sqrt{1 - b^2}$. Let $$\hat{c} = \frac{\sqrt{1-b^2}}{|\vec{c}|}\vec{c}$$ These are uniformly distributed on the sphere of radius $\sqrt{1-b^2}$ in $\mathbb{R}^n$.
  • Use these as the coefficients for the $\vec{u}_i$: the vector $\sum_i^n \hat{c}_i \vec{u}_i$ is now a random point on the unit sphere in the hyperplane $\langle \vec{a}, \vec{x} \rangle = 0$.
  • Translate by $b \vec{a}$ to be on the hyperplane $\langle \vec{a}, \vec{x} \rangle = b$.

Not sure if this is the most efficient algorithm, but it should get the job done. I will hopefully update with some working python code later, using the numpy library to do the random draws.


Python code with explanation

Note that this may not be the most efficient code, but it should get the job done.

n = 5
num_samples = 500
b = 0.7

Colab link to uncommented code. What follows is a line by line explanation of the code:

I will now generate a vector $\vec{a} \in \mathbb{R}^n$ of norm 1:

a = np.random.multivariate_normal(np.zeros(n), np.eye(n), size=1, check_valid='warn', tol=1e-8).reshape(n)
a = a/np.linalg.norm(a)

To project $\vec{x}$ onto the hyperplane $\vec{a}^\perp$, we map $\vec{x} \mapsto \vec{x} - (\vec{x} \cdot \vec{a})\vec{a}$. The matrix of this map is $I - \vec{a} \vec{a}^\top$. This has eigenvalues of $0$ and $1$. $\vec{a}$ is the eigenvector with eigenvalue $0$ and every vector in $\vec{a}^\perp$ is an eigenvector with eigenvalue $1$. The numpy function np.linalg.eigh takes a matrix and returns a tuple. The index 1 element of this tuple is a matrix whose columns are orthonormal eigenvectors, corresponding to increasing eigenvalues. So we want the last $n-1$ of these columns.

U = np.linalg.eigh(np.eye(n) - np.dot(a.reshape(-1,1), a.reshape(1,-1)))[1][:,1:]

I now generate num_sample samples uniformly at random from the unit sphere in $\mathbb{R}^{n-1}$:

c = np.random.multivariate_normal(np.zeros(n-1), np.eye(n-1), size=num_samples, check_valid='warn', tol=1e-8)
c = c/np.linalg.norm(c, axis = 1, keepdims = True)

I use these as the coordinates with the basis $U$ of the subspace $\vec{a}^\perp$:

s = np.dot(U,np.transpose(c)) 

Finally I scale and translate as indicated in my sketch above:

s = s*np.sqrt(1 - b**2)
s = s + np.repeat(b*a.reshape(-1,1), num_samples, axis = 1) 

You can check that these are all on the unit sphere:

np.linalg.norm(s, axis = 0, keepdims = True)

The above should be an array of ones.

You can also check that the dot product of all of these with $\vec{a}$ are equal to $b$:

np.dot(a,s)

The above should be an array whose entries are all $b$.


Visualization

Setting $n = 3$ in the code above gives us the case of intersection a plane with the unit sphere in $3D$ space, which is something we can visualize:

visualization of n=3 case

$\endgroup$
4
  • $\begingroup$ I should have just used np.linalg.eigh instead of np.linalg.svd since $I - \vec{a}\vec{a}^\top$ is symmetric but my head has been in svd land recently. Only slight difference in the code is that eigh has eigenvalues increasing while svd has singular values decreasing. So you would need the last $n-1$ columns instead of the first $n-1$ columns. $\endgroup$ Oct 17, 2023 at 18:36
  • 1
    $\begingroup$ Thanks! Can you share your visualization code as well? $\endgroup$ Oct 17, 2023 at 18:56
  • $\begingroup$ @YaroslavBulatov I included it at the end of the colab link. Hope you don't mind that I sent you a linkedin message as well. I am interested in transitioning to a career in machine learning. $\endgroup$ Oct 17, 2023 at 18:58
  • $\begingroup$ @YaroslavBulatov I also just updated the code to use np.linalg.eigh which should perform better numerically. I also think the intuition is a bit clearer this way. $\endgroup$ Oct 17, 2023 at 19:01

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.