File Exchange

## Minimum Volume Enclosing Ellipsoid

version 1.2.0.0 (2.05 KB) by
Computes the minimum-volume covering ellipoid that encloses N points in a D-dimensional space.

24 Downloads

Updated 20 Jan 2009

View License

[A , c] = MinVolEllipse(P, tolerance)

Finds the minimum volume enclosing ellipsoid (MVEE) of a set of data points stored in matrix P. The following optimization problem is solved:

minimize log(det(A))
s.t. (P_i - c)'*A*(P_i - c)<= 1

in variables A and c, where P_i is the i-th column of the matrix P.
The solver is based on Khachiyan Algorithm, and the final solution is different from the optimal value by the pre-specified amount of 'tolerance'.
---------------------------
Outputs:

c : D-dimensional vector containing the center of the ellipsoid.

A : This matrix contains all the information regarding the shape of the ellipsoid. To get the radii and orientation of the ellipsoid take the Singular Value Decomposition ( svd function in matlab) of the output matrix A:

[U Q V] = svd(A);

the radii are given by:

r1 = 1/sqrt(Q(1,1));
r2 = 1/sqrt(Q(2,2));
...
rD = 1/sqrt(Q(D,D));

and matrix V is the rotation matrix that gives you the orientation of the ellipsoid.

For plotting in 2D or 3D, use MinVolEllipse_plot.m (see the link bellow)

### Cite As

Nima Moshtagh (2021). Minimum Volume Enclosing Ellipsoid (https://www.mathworks.com/matlabcentral/fileexchange/9542-minimum-volume-enclosing-ellipsoid), MATLAB Central File Exchange. Retrieved .

### Comments and Ratings (50)

KUNAL verma

code to get the figure of plot

KUNAL verma

How to plot the MVEE
somebody please share plotting code
its bothering me since 1 month

Jim Lux

Quite useful - I'm finding bounding ellipses for 30 points, and it's plenty fast for that. While I'm sure there are performance optimizations available, I like how this is short, simple, and clear.

(and, it wasn't too tough to turn it into Numpy code, either)

quan yuan

Muhammad Safdar

baopeng liao

thanks, sir. Your code solves a problem that has been bothering me for a long time.

Subhash Nerella

How to get the coordinates of major and minor axis extremeties from the equation obtained?

faride farahnak

thanks, sir. How to get a(semi-major axis) and b (semi-minor axis) for computing the area?

and how to compute the angle between the ellipse and horizon line?

Jeanne Stransky

Wow !! Mind-blowing piece of code ! Thank you so much and also for the explanatory paper !

Nima Moshtagh

I don't check this page frequently.. please send your questions to nima.moshtagh@gmail.com

Mohammad Bhat

How to extract the major and minor axis form (x-c)' * A * (x-c) = 1....so that we will able to draw concentric ellipse within this minimum ellipsoid

Mohammad Bhat

How to draw concentric ellipses...

Mohammad Bhat

Sir , I have to plot minimum ellipsoid , i am also getting same input...my question with some contour coordinates, how to display or plot the calculated minimum eellipsoid , the code given below by Peter Lawrence shows error at U=[cos(thu) -sin(thu);sin(thu) cos(thu)]*[l1*cos(tq);l2*sin(tq)];
Error using *
Inner matrix dimensions must agree.

please help...

BlueEyes

As per script example
--------------
>> format long g
>> P = rand(5,100);
>> [A, c] = MinVolEllipse(P, .01)
A =

Columns 1 through 4:

1.68008932443583 0.451549485542746 0.377295694550125 -0.166024394814601
0.451549485542746 1.88428903019372 0.209271527025232 0.299874743747665
0.377295694550125 0.209271527025232 1.74423212308963 -0.326812488325158
-0.166024394814601 0.299874743747665 -0.326812488325158 1.57099865112731
0.122821795364375 0.0527456078771554 -0.130956861776366 0.0886013717557122

Column 5:

0.122821795364374
0.0527456078771553
-0.130956861776366
0.0886013717557123
1.82006599976003

c =

0.532885969400172
0.543678899108111
0.473317546710603
0.487760240411482
0.423440854549754

>>

Cheers

Mohammad Bhat

How to run it

jianzhong Ding

shuming cheng

Hi Nima,
Thanks for the useful code.

I have a problem when using this code. All of my 3-dimension data is inside the unit Ball, but the output gives a "large" ellipsoid with semi-axes larger than 1. Could you help me find out how can I obtain a more precise result?

Peter Lawrence

I forgot to add to my original comment that, upon calculating [C, m] = MinVolEllipse(data,tol), you need to first invert the matrix. So please use
C=inv(C);
prior to the eigenvalue decomposition.

Peter Lawrence

Works very well for low dimensional data.
@Ali, if you want to plot the ellipse obtained by [C, m] = MinVolEllipse(data,tol), use the following:
tq=linspace(-pi,pi,M); % for display purposes, M is the number of points on the ellipse
[Ve,De]=eig(C);
De=sqrt(diag(De));
[l1,Ie] = max(De);
veig=Ve(:,Ie);
thu=atan2(veig(2),veig(1));
l2=De(setdiff([1 2],Ie));
U=[cos(thu) -sin(thu);sin(thu) cos(thu)]*[l1*cos(tq);l2*sin(tq)];
plot(U(1,:)+m(1),U(2,:)+m(2))

Vincent Phan

Ali Braytee

Ali Braytee

Hi Nima,

How to get the vertices or edge point of the ellipsoid??

`bg`

Hi Nima,

There are a few optimization for your code which could make it much faster and more memory efficient.

Replace X = Q * diag(u) * Q';
by X = bsxfun(@times, Q, u') * Q';

Replace M = diag(Q' * inv(X) * Q);
by M = sum((X\Q).*Q,1);

And finally use:

c = P*u;
A = (1/d) * inv(bsxfun(@times, P, u')*P' - c*c');

Best,
`bg`

Uri Cohen

Thanks for the great code!

I get the following warnings on numerical stability when calculating bounding ellipsoid of ~50 points in ~1000 dimensional space:
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate.
RCOND = 2.952356e-23.
> In MinVolEllipse at 67
In calc_ellipsoid_properties at 33
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate.
RCOND = 1.659718e-21.
> In MinVolEllipse at 89
In calc_ellipsoid_properties at 33

Is the result valid? Is there a way to avoid this?

Nima

The paper on MVEE can be downloaded here: http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.116.7691

Muhammad Anoosh

Thank you for your code. I have two questions:
Question 1 is related to your paper on MVEE and question 2 is regarding the link you mentioned in the comments and ratings of this code.
Q.1 You have mentioned in your paper that the volume of Ellipse is given by v0/(det(sqrt(E)). Can you please elaborate how have you obtained this volume (any reference)

Q.2 You have mentioned a link to the (unpublished) paper describing the algorithm and math behind the code:

The above link is not working can you please give the latest link.

Stefan

Thank you for this code, helps a lot! I am also looking for a method of finding the minimum volume enclosing spheroid, i.e. an ellipsoid where two of the three semi-axes are identical. Do you know of any algorithms for this problem?

David

Hmmm....
det(A) = pi^2 / area^2,

so it seems the optimization for minimum area is to maximize, rather than minimize, log(det(A)).

Also, for small d, it may be easier to obtain semi major/minor lengths (a, b) by finding the eigenvalues/vectors of A (eigenvalues are 1/a^2, 1/b^2, etc, rotation matrix columns are the eigenvectors).

I'm finding this algorithm gives me det(A) < 0 (a hyperbola) for some point sets. Also, A should be invariant under translation of all the input points (which should just result in a translation of c), but I find it is quite sensitive.

Anton Semechko

sweet, thanks!

Pengju

Javier

Hi,

thanks for sharing the code it's really good. One question, please could you explain me which it's the differente between this method and least squares?

Thank you very much

banjo

thanks. This is good.

Sepehr Farhand

Glib

Nima Moshtagh

Here is a link to the (unpublished) paper describing the algorithm and math behind the code:

https://sites.google.com/site/nimamoshtagh/software/MVEE.pdf?attredirects=0

Reza Jafari

I would like to thank you the author for the useful code which he developed. It solved many of my problems. The only problem which I faced was due to the singularity of Matrix X in line 76 and singularity of Matrix P * U * P' - (P * u)*(P*u)' in line 96. Do you have any suggestion to avoid singularity for these two matrices?

Appreciated.

Regards,
Reza

Rashed

Hi,

I have a problem in 3D. Say, I have an 3D array of P(5,5,5) and how can I make this to work? Please suggest. Or, do I need to produce a matrix P containing the coordinates of the region?

Rashed

Thanks a lot. I was looking for something like this.

Nima Moshtagh

Here is a paper describing the algorithm and math behind the code: http://www.seas.upenn.edu/~nima/papers/Mim_vol_ellipse.pdf

Avinash Uppuluri

Can you also provide a good journal reference for this algorithm. Thanks!

Avinash Uppuluri

Hello Nima Moshtagh, The code posted was very useful. Is there a way to add a weight to the points, giving the importance of each point to be included within the ellipse. (Instead of using the tolerance which excludes points through the periphery). With the idea that points that are of less significance need not be included in the ellipse. Where can I include such a weighing parameter. I can probably do it in the pre-processing stage but is there a way to include it within your code. Thanks!

Lien Kuqn Hua

THank you

Santanu Ghorai

Very good program. It has helped me a lot to do with the matlab. Thanks to Nima Moshtagh for submitting the code.

nima moshtagh

Dear Ondrej Danko,
Thank you for bringing this point to my attention. The outcome that a number of points fall outside of the ellipsoid is expected, because we are truncating the algorithm before reaching its optimal value. The amount of error that we will get by doing this, is of the order of parameter 'tolerance', which represents the amount of error you can tolerate in the final result. For example for t=0.001 the average distance of an outside point is of the order of 10^-3.

Ondrej Danko

"Computes the minimum-volume covering ellipoid that encloses N points in a D-dimensional space"...

I would expect that resulting ellipsoid COVERS all points, I mean no points are outside of the resulting ellipsoid...

I had ran a small test if it is so... and surprisingly not all points are necessarily inside.

P = rand(3,N);
[A , C] = MinVolEllipse(P, t);
for j = 1:N
if ((P(:,j)-C)'*A*(P(:,j)-C) > 1)
disp('Point is outside...');
end
end

I am missing something, or it is expected? If expected, any suggestions to modify the algorithm to output ellipsoid strictly covering all points?

Iram Weinstein

I find that the run time can be radically improved by recognizing that the only points that determine the bounding ellipsoid are on the convex hull, as used by John D'Errico in his minboundsphere. Given points P, the following code finds the reduced set of points lying on the hull.

K=convhulln(P.');
K=unique(K(:));
Q=P(:,K);

for P=randn(3,500)
my computer runs MinVolEllipse(P, .001) in 12.54 seconds

and finds Q followed by MinVolEllipse(Q, .001) in 0.08 seconds

M Doug k Vandenberg

This is a very good program for doing minimum volume ellipsoids. it has been very helpful and the author has as well.

Carl Miller

Excellent stuff! This was exactly what I was looking for and works quite well. Although the math behind it all was a bit beyond me, the author was more than helpful in helping me make the most of what the code has to offer, thanks again Nima!

nima moshtagh

For tolerance=0.01 in 2D here is how long it takes (approximatly) to find the covering ellipse on a 2.4GHz computer:

#point time(sec)
------ --------
100 0.03
500 1.24
1000 5.20
1500 13.26
2000 25.27

John D'Errico

This could be useful IF you have small sets of data.

I tested this for several problems, in 1, 2 and 5 dimensions. All seemed to work. The only downside is data of even moderate size. The author gives an example with 100 points, solved fairly quickly. But expand that to only 500 points (in 2-d) and it takes 140 seconds to run on my computer.

##### MATLAB Release Compatibility
Created with R13
Compatible with any release
##### Platform Compatibility
Windows macOS Linux

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!