def mvee(points, tol=0.001):
    """
    Returns the minimum-volume enclosing ellipse (MVEE)
    of a set of points, using the Khachiyan algorithm.
    """

    # This function is a port of the matlab function by
    # Nima Moshtagh found here:
    # https://www.mathworks.com/matlabcentral/fileexchange/9542-minimum-volume-enclosing-ellipsoid
    # with accompanying writup here:
    # https://www.researchgate.net/profile/Nima_Moshtagh/publication/254980367_MINIMUM_VOLUME_ENCLOSING_ELLIPSOIDS/links/54aab5260cf25c4c472f487a.pdf

    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.outer(c,c) )
    return A, np.transpose(c)
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 . outer ( c , c ) ) return A , np . transpose ( c )
Returns the minimum-volume enclosing ellipse (MVEE)
    of a set of points, using the Khachiyan algorithm.