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