Drawing uniform random numbers from within an ellipsoid

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. Thank you so much for the article Matthew. Really helped me a lot 🙂

   Reply

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

   Reply

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
   

   Reply

4. 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

   Reply

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

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

      Reply

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

   Reply