Drawing samples from a multivariate normal distribution subject to quadratic constraints I would like to efficiently draw samples $x \in \mathbb{R}^d$ from $\mathcal{N}(\mu, \Sigma)$ subject to the constraint that $||x||_2 = 1$.
 A: The formal resolution of this problem first requires a proper definition of a 

"$\mathcal{N}_d(μ,Σ)$ distribution subject to the constraint that
  $||x||^2=1$"

The natural way is to define the distribution of $X\sim\mathcal{N}_d(μ,Σ)$ conditional on $||X||=\varrho$. And to apply this conditional to the case $\varrho=1$. If we use polar coordinates, 
$$\eqalign{
x_1&=\varrho\cos(\theta_1)\qquad&\theta_1\in[0,\pi]\\
x_2&=\varrho\sin(\theta_1)\cos(\theta_2)\qquad&\theta_2\in[0,\pi]\\
&\vdots\\
x_{d-1}&=\varrho \left( \prod_{i=1}^{d-2}\sin(\theta_i) \right) \cos(\theta_{d-1})\qquad&\theta_{d-1}\in[0,2\pi]\\
x_{d}&=\varrho\prod_{i=1}^{d-1}\sin(\theta_i)
}$$
the Jacobian of the transform is
$$\varrho^{d-1}\prod_{i=1}^{d-2}\sin(\theta_i)^{d-1-i}$$
Therefore the conditional density of the distribution of $\mathbf{\theta}=(\theta_1,\ldots,\theta_{d-1})$ given $\varrho$ is
$$
f(\mathbf{\theta}|\varrho) \propto \exp\frac{-1}{2}\left\{(x(\theta,\varrho)-\mu)^T\Sigma^{-1}(x(\theta,\varrho)-\mu) \right\} \prod_{i=1}^{d-2}\sin(\theta_i)^{d-1-i}$$

Conclusion: This density differs from simply applying the Normal density to a point on the unit sphere because of the Jacobian.

The second step is to consider the target density
$$
f(\mathbf{\theta}|\varrho=1) \propto \exp\frac{-1}{2}\left\{(x(\theta,1)-\mu)^T\Sigma^{-1}(x(\theta,1)-\mu) \right\} \prod_{i=1}^{d-2}\sin(\theta_i)^{d-1-i}$$
and design a Markov chain Monte Carlo algorithm to explore the parameter space $[0,\pi]^{d-2}\times[0,2\pi]$. My first attempt would be at a Gibbs sampler, initialised at the point on the sphere closest to $\mu$, that is, $\mu/||\mu||$, and proceeding one angle at a time in a Metropolis-within-Gibbs manner:


*

*Generate $\theta_1^{(t+1)}\sim\mathcal{U}([\theta_1^{(t)}-\delta_1,\theta_1^{(t)}+\delta_1])$ (where sums are computed modulo $\pi$) and accept this new value with probability
$$\dfrac{f(\theta_1^{(t+1)},\theta_2^{(t)},...|\varrho=1)}{f(\theta_1^{(t)},\theta_2^{(t)},...|\varrho=1)}\wedge 1$$
else $\theta_1^{(t+1)}=\theta_1^{(t)}$

*Generate $\theta_2^{(t+1)}\sim\mathcal{U}([\theta_2^{(t)}-\delta_2,\theta_2^{(t)}+\delta_2])$ (where sums are computed modulo $\pi$) and accept this new value with probability
$$\dfrac{f(\theta_1^{(t+1)},\theta_2^{(t+1)},\theta_3^{(t)},...|\varrho=1)}{f(\theta_1^{(t+1)},\theta_2^{(t)},\theta_3^{(t)},...|\varrho=1)}\wedge 1$$
else $\theta_2^{(t+1)}=\theta_2^{(t)}$

*$\ldots$

*Generate $\theta_{d-1}^{(t+1)}\sim\mathcal{U}([\theta_{d-1}^{(t)}-\delta_{d-1},\theta_{d-1}^{(t)}+\delta_{d-1}])$ (where sums are computed modulo $2\pi$) and accept this new value with probability
$$\dfrac{f(\theta_1^{(t+1)},\theta_2^{(t+1)},...,\theta_{d-1}^{(t+1)}|\varrho=1)}{f(\theta_1^{(t+1)},\theta_2^{(t+1)},...,\theta_{d-1}^{(t)}|
\varrho=1)}\wedge 1$$
else $\theta_{d-1}^{(t+1)}=\theta_{d-1}^{(t)}$


The scales $\delta_1$, $\delta_2$, $\ldots$, $\delta_{d-1}$ can be scaled against the acceptance rates of the steps, towards an ideal goal of $50\%$.
Here is an R code to illustrate the above, with default values for $\mu$ and $\Sigma$:
library(mvtnorm)
d=4
target=function(the,mu=1:d,sigma=diag(1/(1:d))){
 carte=cos(the[1])
 for (i in 2:(d-1))
  carte=c(carte,prod(sin(the[1:(i-1)]))*cos(the[i]))
 carte=c(carte,prod(sin(the[1:(d-1)])))
 prod(sin(the)^((d-2):0))*dmvnorm(carte,mean=mu,sigma=sigma)}
#Gibbs
T=1e4
#starting point
mu=(1:d)
mup=mu/sqrt(sum(mu^2))
mut=acos(mup[1])
for (i in 2:(d-1))
  mut=c(mut,acos(mup[i]/prod(sin(mut))))
thes=matrix(mut,nrow=T,ncol=d-1,byrow=TRUE)
delta=rep(pi/2,d-1)     #scale
past=target(thes[1,])   #current target
for (t in 2:T){
 thes[t,]=thes[t-1,]
 for (j in 1:(d-1)){
   prop=thes[t,]
   prop[j]=prop[j]+runif(1,-delta[j],delta[j])
   prop[j]=prop[j]%%(2*pi-(j<d-1)*pi)
   prof=target(prop)
   if (runif(1)<prof/past){
     past=prof;thes[t,]=prop}
   }
}

A: $||x||_2^2=1$ is not strictly possible since $x$ is a (continuous) random variable. If you would like it to have a variance of 1, i.e. $E[(x-\mu)^2]\tilde{=} \frac{1}{n}\sum (x-\mu)^2=\frac{1}{n} ||x-n||_2^2=\frac{1}{n}$ (where the tilde means we estimate the variance), then you would need to require its variance to be $\frac{1}{n}$. However, this demand may conflict with $\Sigma$. That is, to get samples with this variance you need the diagonal of $\Sigma$ to be equal to $\frac{1}{n}$.
To sample form this distribution in general, you can generate i.i.d standard normals, and then multiply by $\Sigma^{0.5}$, the square root of $\Sigma$, and then add the means $\mu$. 
