# 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. benlambert18785 says:

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

# 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

1. Felix says:

It seems similar to what is done here, even though I was hoping for something even simpler:

https://ieeexplore.ieee.org/document/758215

5. G Prathap says:

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;
}