MATLAB code exists to find the so-called "minimum volume enclosing ellipsoid" (e.g. here, also here). I'll paste the relevant part for convenience:
function [A , c] = MinVolEllipse(P, tolerance)
[d N] = size(P);
Q = zeros(d+1,N);
Q(1:d,:) = P(1:d,1:N);
Q(d+1,:) = ones(1,N);
count = 1;
err = 1;
u = (1/N) * ones(N,1);
while err > tolerance,
X = Q * diag(u) * Q';
M = diag(Q' * inv(X) * Q);
[maximum j] = max(M);
step_size = (maximum - d -1)/((d+1)*(maximum-1));
new_u = (1 - step_size)*u ;
new_u(j) = new_u(j) + step_size;
count = count + 1;
err = norm(new_u - u);
u = new_u;
end
U = diag(u);
A = (1/d) * inv(P * U * P' - (P * u)*(P*u)' );
c = P * u;
Here is some MATLAB test code:
points = [[ 0.53135758, -0.25818091, -0.32382715],
[ 0.58368177, -0.3286576, -0.23854156,],
[ 0.18741533, 0.03066228, -0.94294771],
[ 0.65685862, -0.09220681, -0.60347573],
[ 0.63137604, -0.22978685, -0.27479238],
[ 0.59683195, -0.15111101, -0.40536606],
[ 0.68646128, 0.0046802, -0.68407367],
[ 0.62311759, 0.0101013, -0.75863324]];
[A centroid] = minVolEllipse(points',0.001);
A
[~, D, V] = svd(A);
rx = 1/sqrt(D(1,1));
ry = 1/sqrt(D(2,2));
rz = 1/sqrt(D(3,3));
[u v] = meshgrid(linspace(0,2*pi,20),linspace(-pi/2,pi/2,10));
x = rx*cos(u').*cos(v');
y = ry*sin(u').*cos(v');
z = rz*sin(v');
for idx = 1:20,
for idy = 1:10,
point = [x(idx,idy) y(idx,idy) z(idx,idy)]';
P = V * point;
x(idx,idy) = P(1)+centroid(1);
y(idx,idy) = P(2)+centroid(2);
z(idx,idy) = P(3)+centroid(3);
end
end
figure
plot3(points(:,1),points(:,2),points(:,3),'.');
hold on;
mesh(x,y,z);
axis square;
alpha 0;
which will produce the the covariance matrix:
A =
47.3693 -116.0758 -79.1861
-116.0758 458.0874 280.0656
-79.1861 280.0656 179.3886
Now, here is my attempt at port this code to Python (2.7):
from __future__ import division
import numpy as np
import numpy.linalg as la
def mvee(points,tol=0.001):
N, d = points.shape
Q = np.zeros([N,d+1])
Q[:,0:d] = points[0:N,0:d]
Q[:,d] = np.ones([1,N])
Q = np.transpose(Q)
points = np.transpose(points)
count = 1
err = 1
u = (1/N) * np.ones(shape = (N,))
while err > tol:
X = np.dot(np.dot(Q, np.diag(u)), np.transpose(Q))
M = np.diag( np.dot(np.dot(np.transpose(Q), la.inv(X)),Q))
jdx = np.argmax(M)
step_size = (M[jdx] - d - 1)/((d+1)*(M[jdx] - 1))
new_u = (1 - step_size)*u
new_u[jdx] = new_u[jdx] + step_size
count = count + 1
err = la.norm(new_u - u)
u = new_u
U = np.diag(u)
c = np.dot(points,u)
A = (1/d) * la.inv(np.dot(np.dot(points,U), np.transpose(points)) - np.dot(c,np.transpose(c)) )
return A, np.transpose(c)
The corresponding test code:
from mpl_toolkits.mplot3d import Axes3D
from mpl_toolkits.mplot3d.art3d import Poly3DCollection
import matplotlib.pyplot as plt
from scipy.spatial import Delaunay
#some random points
points = np.array([[ 0.53135758, -0.25818091, -0.32382715],
[ 0.58368177, -0.3286576, -0.23854156,],
[ 0.18741533, 0.03066228, -0.94294771],
[ 0.65685862, -0.09220681, -0.60347573],
[ 0.63137604, -0.22978685, -0.27479238],
[ 0.59683195, -0.15111101, -0.40536606],
[ 0.68646128, 0.0046802, -0.68407367],
[ 0.62311759, 0.0101013, -0.75863324]])
# compute mvee
A, centroid = mvee(points)
print A
# point it and some other stuff
U, D, V = la.svd(A)
rx, ry, rz = [1/np.sqrt(d) for d in D]
u, v = np.mgrid[0:2*np.pi:20j,-np.pi/2:np.pi/2:10j]
x=rx*np.cos(u)*np.cos(v)
y=ry*np.sin(u)*np.cos(v)
z=rz*np.sin(v)
for idx in xrange(x.shape[0]):
for idy in xrange(y.shape[1]):
x[idx,idy],y[idx,idy],z[idx,idy] = np.dot(np.transpose(V),np.array([x[idx,idy],y[idx,idy],z[idx,idy]])) + centroid
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')
ax.scatter(points[:,0],points[:,1],points[:,2])
ax.plot_surface(x, y, z, cstride = 1, rstride = 1, alpha = 0.1)
plt.show()
produces this:
[[ 0.84650504 -1.40006147 0.39857055]
[-1.40006147 2.60678264 -1.52583781]
[ 0.39857055 -1.52583781 1.04581752]]
Clearly different. What gives?
