# Drawing uniform random numbers from within an ellipsoid

Here’s a Matlab function I wrote to draw a set of random numbers from a uniform distribution within an n-dimensional ellipsoid, where the size and orientation of the ellipsoid is defined by a covariance matrix.

function pnts = draw_from_ellipsoid( covmat, cent, npts )

% function pnts = draw_from_ellipsoid( covmat, cent, npts )
%
% This function draws points uniformly from an n-dimensional ellipsoid
% with edges and orientation defined by the the covariance matrix covmat.
% The number of points produced in the n-dimensional space is given by
% npts, with an output array of npts by ndims. The centre of each dimension
% is given in the vector cent.

% get number of dimensions of the covariance matrix
ndims = length(covmat);

% calculate eigenvalues and vectors of the covariance matrix
[v, e] = eig(covmat);

% check size of cent and transpose if necessary
sc = size(cent);
if sc(1) > 1
cent = cent';
end

% generate radii of hyperspheres
rs = rand(npts,1);

% generate points
pt = randn(npts,ndims);

% get scalings for each point onto the surface of a unit hypersphere
fac = sum(pt(:,:)'.^2);

% calculate scaling for each point to be within the unit hypersphere
% with radii rs
fac = (rs.^(1/ndims)) ./ sqrt(fac');

pnts = zeros(npts,ndims);

% scale points to the ellipsoid using the eigenvalues and rotate with
% the eigenvectors and add centroid
d = sqrt(diag(e));
for i=1:npts
% scale points to a uniform distribution within unit hypersphere
pnts(i,:) = fac(i)*pt(i,:);

% scale and rotate to ellipsoid
pnts(i,:) = (pnts(i,:) .* d' * v') + cent;
end

end


## 6 Replies to “Drawing uniform random numbers from within an ellipsoid”

1. Krishna says:

Thank you so much for the article Matthew. Really helped me a lot 🙂

2. John A. Major says:

Thanks – I needed that and translated it into R. Will send you code if you want.

3. Thanks very much for this. In case anyone wants the code for Python, here you go:

def draw_from_ellipsoid( covmat, cent, npts):
# random uniform points within ellipsoid as per: https://ma.ttpitk.in/blog/?p=368
ndims = covmat.shape[0]

# calculate eigenvalues (e) and eigenvectors (v)
eigenValues, eigenVectors = np.linalg.eig(covmat)
idx = (-eigenValues).argsort()[::-1][:ndims]
e = eigenValues[idx]
v = eigenVectors[:,idx]
e = np.diag(e)

# generate radii of hyperspheres
rs = np.random.uniform(0,1,npts)

# generate points
pt = np.random.normal(0,1,[npts,ndims]);

# get scalings for each point onto the surface of a unit hypersphere
fac = np.sum(pt**2,axis=1)

# calculate scaling for each point to be within the unit hypersphere
# with radii rs
fac = (rs**(1.0/ndims)) / np.sqrt(fac)

pnts = np.zeros((npts,ndims));

# scale points to the ellipsoid using the eigenvalues and rotate with
# the eigenvectors and add centroid
d = np.sqrt(np.diag(e))
d.shape = (ndims,1)

for i in range(0,npts):
# scale points to a uniform distribution within unit hypersphere
pnts[i,:] = fac[i]*pt[i,:]
pnts[i,:] = np.dot(np.multiply(pnts[i,:],np.transpose(d)),np.transpose(v)) + cent

return pnts

4. Felix says:

Hi, thanks a lot! Can you point me to the underlying theory?

I get the intuition:

You sample points from a normal distribution and normalize them to be on the unit sphere. That way, you get points with a uniform distribution in their orientation. Then you combine the uniform orientations with uniform radii (is that results still uniform?). Where is the exponent in r^*(1/ndim) coming from?
Finally, you rotate and rescale.

Looking forward to some hint,

Felix

5. Someone needs this in c++, here you go..

Eigen::MatrixXf generate_samples_from_ellipsoid(Eigen::MatrixXf covmat, Eigen::VectorXf cent, int npts){
int ndims = covmat.rows();
Eigen::EigenSolver eigensolver;
eigensolver.compute(covmat);
Eigen::VectorXf eigen_values = eigensolver.eigenvalues().real();
Eigen::MatrixXf eigen_vectors = eigensolver.eigenvectors().real();
std::vector<std::tuple> eigen_vectors_and_values;

for(int i=0; i<eigen_values.size(); i++){
std::tuple vec_and_val(eigen_values[i], eigen_vectors.row(i));
eigen_vectors_and_values.push_back(vec_and_val);
}
std::sort(eigen_vectors_and_values.begin(), eigen_vectors_and_values.end(),
[&](const std::tuple& a, const std::tuple& b) -> bool{
return std::get(a) <= std::get(b);
});
int index = 0;
for(auto const vect : eigen_vectors_and_values){
eigen_values(index) = std::get(vect);
eigen_vectors.row(index) = std::get(vect);
index++;
}

Eigen::MatrixXf eigen_values_as_matrix = eigen_values.asDiagonal();

std::random_device rd{};
std::mt19937 gen{rd()};
std::uniform_real_distribution dis(0, 1);
std::normal_distribution normal_dis{0.0f, 1.0f};

Eigen::MatrixXf pt = Eigen::MatrixXf::Zero(npts, ndims).unaryExpr([&](float dummy){return (float)normal_dis(gen);});
Eigen::VectorXf rs = Eigen::VectorXf::Zero(npts).unaryExpr([&](float dummy){return dis(gen);});
Eigen::VectorXf fac = pt.array().pow(2).rowwise().sum();
Eigen::VectorXf fac_sqrt = fac.array().sqrt();
Eigen::VectorXf rs_pow = rs.array().pow(1.0/ndims);
fac = rs_pow.array()/fac_sqrt.array();
Eigen::MatrixXf points = Eigen::MatrixXf::Zero(npts, ndims);
Eigen::VectorXf d = eigen_values_as_matrix.diagonal().array().sqrt();
for(auto i(0); i<npts; i++){
points.row(i) = fac(i)*pt.row(i).array();
Eigen::MatrixXf fff = (points.row(i).array()*d.transpose().array());
Eigen::VectorXf bn = eigen_vectors*fff.transpose();
points.row(i) = bn.array() + cent.array();
}
std::cout << "points: " << points << std::endl;
return points;
}