Section 3.4.1: Numerical Distributions

Exercise 28. [HM35] (R. P. Brent.) Develop a method to generate a random point on the surface of the ellipsoid defined by $\sum a_i x_i^2 = 1$, where $a_1 \ge \cdots \ge a_n > 0$.

Verified: yes
Solve time: 1m43s

---

Setup

Let

$$ \qquad a_1\ge \cdots \ge a_n>0. $$

The problem is to generate a random point on the surface of the ellipsoid $E$, with respect to the natural surface-area measure.

Write

$$$$

so that $E={x:x^{T}Ax=1}$.

The task is to construct an algorithm whose output distribution on $E$ is proportional to surface area.

Solution

Let $Y=(Y_1,\ldots,Y_n)$ be a vector of independent normal deviates with mean $0$ and variance $1$. By rotational symmetry,

$$ U=\frac{Y}{\sqrt{Y^{T}Y}} $$

is uniformly distributed on the unit sphere

$$ S^{n-1}={u:u^{T}u=1}. $$

Define the linear transformation

$$ T=A^{-1/2} =\operatorname{diag}!\left( \frac1{\sqrt{a_1}},\ldots, \frac1{\sqrt{a_n}} \right). $$

Then $T(S^{n-1})=E$, since

$$ (Tu)^TA(Tu) =u^TT^TATu =u^TIu =1. $$

A uniformly distributed point $U$ on the sphere does not map under $T$ to a point uniformly distributed with respect to surface area on $E$. The surface element changes by a factor that must be determined.

Let $u\in S^{n-1}$, and let $d\sigma(u)$ denote spherical surface measure. If $dS(x)$ denotes surface measure on $E$, the area formula for a linear map gives

$$ dS(Tu) = |\det T|, |T^{-T}u|, d\sigma(u). $$

Since

$$ T^{-T}=A^{1/2}, $$

the density of the image of uniform spherical measure relative to surface area on $E$ is proportional to

$$ \frac1{|A^{1/2}u|}. $$

Because $u^{T}u=1$,

$$ |A^{1/2}u| = \left(\sum_{i=1}^{n}a_i u_i^2\right)^{1/2}. $$

Hence, if $U$ is uniform on the sphere, the point $X=TU$ has surface-density proportional to

$$ \left(\sum_{i=1}^{n}a_iU_i^2\right)^{-1/2}. $$

To obtain the uniform surface distribution on $E$, compensate for this factor by rejection.

Let

$$ M=\sqrt{a_1}. $$

Since $a_1\ge a_i$,

$$ \left(\sum_{i=1}^{n}a_iU_i^2\right)^{1/2} \le \sqrt{a_1}=M. $$

The following procedure therefore has acceptance probability at most $1$.

Generate $U$ uniformly on $S^{n-1}$. Generate an independent uniform deviate $V$. Accept $U$ if

$$ V\le \frac{\left(\sum_{i=1}^{n}a_iU_i^2\right)^{1/2}} {\sqrt{a_1}}, $$

otherwise generate a new $U$.

After acceptance, output

$$ X=TU = \left( \frac{U_1}{\sqrt{a_1}}, \ldots, \frac{U_n}{\sqrt{a_n}} \right). $$

The accepted spherical points have density proportional to

\left(\sum_{i=1}^{n}a_iU_i^2\right)^{1/2} ] with respect to (d\sigma). Multiplying this by the image-density factor

\left(\sum_{i=1}^{n}a_iU_i^2\right)^{-1/2}

produces a constant density on (E). Therefore the output point (X) is distributed according to surface area on the ellipsoid. Thus the algorithm is: 1. Generate independent standard normal deviates (Y_1,\ldots,Y_n). 2. Set

U_i=\frac{Y_i}{\sqrt{\sum_jY_j^2}}.

3. Generate (V\sim U(0,1)). 4. If

V>

\frac{\left(\sum_i a_iU_i^2\right)^{1/2}}

{\sqrt{a_1}},

$$ return to step 1. 5. Output $$

X_i=\frac{U_i}{\sqrt{a_i}}

\qquad(1\le i\le n).

This generates a random point uniformly distributed over the surface of the ellipsoid. ## Verification The output lies on the ellipsoid because

\sum_{i=1}^{n}a_iX_i^2

\sum_{i=1}^{n}a_i

\left(\frac{U_i}{\sqrt{a_i}}\right)^2

\sum_{i=1}^{n}U_i^2

For a linear map (T), the induced surface-area element satisfies

dS(Tu)

|\det T|,|T^{-T}u|,d\sigma(u).

Substituting (T=A^{-1/2}) yields

dS(Tu)

|\det A^{-1/2}|

\left(\sum_i a_i u_i^2\right)^{1/2}

d\sigma(u).

$$ Hence the image of uniform spherical measure has density proportional to $$

\left(\sum_i a_i u_i^2\right)^{-1/2}

$$ relative to surface area. Rejection with acceptance probability proportional to $$

\left(\sum_i a_i u_i^2\right)^{1/2}

exactly cancels this factor, leaving constant density on the ellipsoid. Therefore the accepted points are uniformly distributed with respect to surface area. This completes the proof. ∎ ## Notes If the requirement were to generate points according to the image of the uniform distribution on the sphere under the map (x=T u), no rejection step would be needed. The rejection step is necessary because the problem asks for uniformity on the ellipsoid surface itself, not on its spherical preimage.