使用Octave，我发现MinVolEllipse中的while循环结束后，
u = 0.0053531 0.2384227 0.2476188 0.0367063 0.0257947 0.2124423 0.0838103 0.1498518
这与
uPython函数找到的结果一致
mvee。在八度音阶上产生更多的调试打印语句
(P*u) = 0.50651 -0.11166 -0.57847
和
(P*u)*(P*u)' = 0.256555 -0.056556 -0.293002 -0.056556 0.012467 0.064590 -0.293002 0.064590 0.334628
但是在Python方面，
c = np.dot(points.T,u)print(c)
产量
[ 0.50651212 -0.11165724 -0.57847018]
和
print(np.dot(c,np.transpose(c)))
产量
0.60364961984 # <-- This should equal (P*u)*(P*u)', a 3x3 matrix.
一旦知道了问题，解决方案就很简单。
(P*u)*(P*u)'可以用以下公式计算：
np.multiply.outer(c,c)
import numpy as npimport numpy.linalg as laimport matplotlib.pyplot as pltfrom mpl_toolkits.mplot3d import Axes3Dpi = np.pisin = np.sincos = np.cosdef mvee(points, tol = 0.001): """ Finds the ellipse equation in "center form" (x-c).T * A * (x-c) = 1 """ N, d = points.shape Q = np.column_stack((points, np.ones(N))).T err = tol+1.0 u = np.ones(N)/N while err > tol: # assert u.sum() == 1 # invariant X = np.dot(np.dot(Q, np.diag(u)), Q.T) M = np.diag(np.dot(np.dot(Q.T, la.inv(X)), Q)) jdx = np.argmax(M) step_size = (M[jdx]-d-1.0)/((d+1)*(M[jdx]-1.0)) new_u = (1-step_size)*u new_u[jdx] += step_size err = la.norm(new_u-u) u = new_u c = np.dot(u,points) A = la.inv(np.dot(np.dot(points.T, np.diag(u)), points) - np.multiply.outer(c,c))/d return A, c#some random pointspoints = 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]])# Singular matrix error!# points = np.eye(3)A, centroid = mvee(points) U, D, V = la.svd(A) rx, ry, rz = 1./np.sqrt(D)u, v = np.mgrid[0:2*pi:20j, -pi/2:pi/2:10j]def ellipse(u,v): x = rx*cos(u)*cos(v) y = ry*sin(u)*cos(v) z = rz*sin(v) return x,y,zE = np.dstack(ellipse(u,v))E = np.dot(E,V) + centroidx, y, z = np.rollaxis(E, axis = -1)fig = plt.figure()ax = fig.add_subplot(111, projection='3d')ax.plot_surface(x, y, z, cstride = 1, rstride = 1, alpha = 0.05)ax.scatter(points[:,0],points[:,1],points[:,2])plt.show()
顺便说一下，此计算使用了大量的矩阵乘法，使用时
np.dot看起来很冗长。如果将NumPy数组转换为NumPy矩阵，则矩阵乘法可以用表示
*。例如，
A = la.inv(np.dot(np.dot(points.T, np.diag(u)), points)- np.dot(c[:, np.newaxis], c[np.newaxis, :]))/d
变成
A = la.inv(points.T*np.diag(u)*points - c.T*c)/d
由于可读性很重要，因此您可能希望使用NumPy矩阵进行主要计算：
def mvee(points, tol = 0.001): """ Find the minimum volume ellipse. Return A, c where the equation for the ellipse given in "center form" is (x-c).T * A * (x-c) = 1 """ points = np.asmatrix(points) N, d = points.shape Q = np.column_stack((points, np.ones(N))).T err = tol+1.0 u = np.ones(N)/N while err > tol: # assert u.sum() == 1 # invariant X = Q * np.diag(u) * Q.T M = np.diag(Q.T * la.inv(X) * Q) jdx = np.argmax(M) step_size = (M[jdx]-d-1.0)/((d+1)*(M[jdx]-1.0)) new_u = (1-step_size)*u new_u[jdx] += step_size err = la.norm(new_u-u) u = new_u c = u*points A = la.inv(points.T*np.diag(u)*points - c.T*c)/d return np.asarray(A), np.squeeze(np.asarray(c))
