I have previously written about how to efficiently generate points uniformly at random on the surface of a sphere. The method uses a mathematical fact from multivariate statistics: If X is drawn from the uncorrelated multivariate normal distribution, then S = X / ||X|| has the uniform distribution on the unit sphere. This is true in any dimension; see the previous article for pictures of a uniform sample on the circle in the plane and on a sphere in 3-D.

It is worth mentioning that this procedure can also be used to generate random points uniformly on any ellipse. This article shows the mathematical algorithm and implements it in SAS. Examples are shown in 2-D (an ellipse in the plane) and in 3-D (an ellipsoid). The same process works for hyperellipsoids in high-dimensional Euclidean space.

From spheres to ellipsoids

Intuitively, you can move from spheres to ellipsoids by changing the covariance matrix that you use for the multivariate normal (MVN) distribution. When the covariance matrix is the identity matrix, the equi-probability contours of the density are spheres. However, for a general positive definite covariance matrix, the equi-probability contours of the density are (hyper)ellipsoids. A particularly simple choice for a covariance matrix is a diagonal matrix. For a diagonal matrix, the axes of the density ellipsoids are orthogonal and are aligned with the standard coordinate axes.

There is a simple geometric relationship between the equi-probability contours for an uncorrelated MVN distribution and for a correlated MVN distribution with covariance matrix Σ. Briefly, if you use the Cholesky factorization to decompose Σ = LL^T, then the matrix L is a linear transformation that maps the spherical contours of the uncorrelated distribution to corresponding contours of the correlated distribution, as shown by the image to the right.

The geometry (and the computation) is simplest when you use a diagonal covariance matrix. If $\Sigma = \mbox{diag}(\sigma_1^2, \sigma_2^2, \ldots, \sigma_d^2)$ is a diagonal covariance matrix, then L = diag(σ₁, σ₂, ..., σ_d) for σ_i > 0. The linear transformation, L, maps the unit sphere at the origin to an ellipsoid with axes that are aligned with the standard coordinate axes and whose radii are σ₁, σ₂, ..., σ_d.

The uniform distribution is preserved under a linear transformation

Of course, this problem requires more than mapping a sphere to an ellipsoid. We want to generate uniform random variates on the ellipsoid. Fortunately, you can prove that the uniform distribution is preserved under a linear transformation. More precisely, suppose that X is a continuous random variable with a uniform distribution on some subset S in d-dimensional Euclidean space. Let Y = μ + L*X be any affine transformation of X, where μ is a d-dimensional vector and L is an invertible d x d matrix. Then Y is uniformly distributed on the image of S, which is the set {y = μ + L*x | x ∈ S}.

This theorem tells you that you can generate random uniform variates on the unit sphere (as shown in the previous article) and then linearly transform those variates to obtain uniform variates on an ellipsoid. The transformed variates will be uniformly distributed on the image of the sphere, which is an ellipsoid.

A SAS program to generate random uniform variates on an ellipse

You can carry out this process in any matrix-vector language, such as R, MATLAB, or SAS IML. The following SAS IML program defines a function that takes two arguments: the number of points (n) to generate, and a vector of d positive values, which represent the elements of a diagonal covariance matrix. It returns n points in d-dimensional space that are on the ellipse {y = diag(σ₁, σ₂, ..., σ_d)*x | x is on the unit sphere and σ_i > 0}.

/* Uniform random points on ellipsoid in d-dimensional space.
   The ellipsoid is the image of the unit sphere under the transformation x -> L*x
   where x is on the unit sphere and L = diag(sigma_1, sigma_2, ..., sigma_d), sigma_i > 0.
*/
proc iml;
/* generate n random uniform points on the ellipse L*x, where x is on unit sphere and L = diag(s) */
start RandEllipse(n, s);
   d = ncol(s);
   mu = j(1, d, 0);
   X = randnormal(n, mu, I(d));       /* X ~ N(0, I) */
   Z = X / sqrt( X[,##] );            /* Z is uniform on circle/sphere */
   L = diag(s);                       /* Cholesky root of covariance matrix Sigma=diag(s##2) */
   return Z * L`;                     /* points are uniform on ellipse/ellipsoid */
finish;
 
call randseed(12345);
/* 2-D example with n=200 points on ellipse in plane */
n = 200;
sigmas = {4  3};
Y = RandEllipse(n, sigmas);
title "200 Random Points on Ellipse";
call scatter(Y[,1], Y[,2]) grid={x y}
             procopt="aspect=0.75" option="transparency=0.5";

A scatter plot of the output shows that the 200 points are distributed uniformly at random on the ellipse with one axis aligned with the vector (4,0) and the other axis aligned with the vector (0,3).

The same function supports generating points on higher-dimensional ellipsoids. For example, to get 1,000 random points on an ellipsoid in 3-D, you can use the following statements.

/* 3-D example with n=1,000 points on ellipsoid */
Y = RandEllipse(1000, {4  3  1});     /* axes aligned with (4,0,0), (0,3,0), and (0,0,1) */
 
/* write points to a data set and visualize */
create RandEllipsoid from Y[c={'x1' 'x2' 'x3'}];
   append from Y;
close;
QUIT;
 
title "1,000 Random Points on an Ellipsoid";
proc sgplot data=RandEllipsoid;
   scatter x=x1 y=x2 / colorresponse=x3 transparency=0.5 markerattrs=(symbol=CircleFilled);
   xaxis grid;
   yaxis grid;
run;

The graph projects the points onto the (x,y) plane. The markers are color-coded according to the z value, which ranges in [-1,1]. Because the points are on an ellipsoid, points near (x,y)=(0,0) have z coordinates near ±1. Points whose values are close to z=0 are located near the ellipse that passes through (±4, 0) and (0, ±3).

Ellipsoids not aligned with the coordinate axis

The same technique works even if you do not use a diagonal covariance matrix. Given any positive definite matrix, you can find the Cholesky decomposition Σ = LL^T and use L to map random variate from the circle onto an ellipse. The following SAS IML function implements this general technique:

 
proc iml;
/* generate n random uniform points on the ellipse L*x, where x is on unit sphere and 
   L satisfies Sigma = L*L` */
start RandEllipseCov(n, Sigma);
   d = ncol(Sigma);
   mu = j(1, d, 0);
   X = randnormal(n, mu, I(d));       /* X ~ N(0, I) */
   Z = X / sqrt( X[,##] );            /* Z is uniform on circle/sphere */
   U = root(Sigma);                   /* Cholesky decomp Sigma=U`*U, U upper triangular */
   L = U`;                            /* Sigma=L*L', L lower triangular */
   return Z * L`;                     /* points are uniform on ellipse/ellipsoid */
finish;
 
call randseed(12345);
/* 2-D example with n=200 points on ellipse in plane */
n = 200;
Sigma = {16  4,
          4  9};
Y = RandEllipseCov(n, Sigma);
 
title "200 Random Points on Tilted Ellipse";
/* slopes of reference lines are the directions of the eigenvectors of Sigma */
call scatter(Y[,1], Y[,2]) grid={x y}
             procopt="aspect=0.75 noautolegend" option="transparency=0.5"
             other="lineparm x=0 y=0 slope=0.4538/ lineattrs=(color=gray); lineparm x=0 y=0 slope=-2.2038/lineattrs=(color=gray);";

The graph overlays two reference lines, which indicate the axes for the ellipse. They are computed by using the eigenvectors of the Σ matrix. The computation is not shown here, but you can read the details in a previous article about ellipses and covariance matrices.

Summary

This article shows how to generate points uniformly at random on an ellipse, on an ellipsoid, or on any hyperellipsoid in higher dimensions. It first generates points uniformly at random on the unit sphere, then uses a linear transformation to map the points onto an ellipsoid. If you want an ellipse whose axes are aligned with the coordinate axes, the linear transformation is given by a simple diagonal matrix with positive values on the diagonal. You can also generate points on ellipsoids whose axes are not aligned with the coordinate axes. If Σ is any positive definite matrix, you can decompose Σ = LL^T. If you generate points on the unit sphere, you can use L to map them onto the ellipse whose axes are aligned with the eigendirections of Σ.