ellipsoid.py 文件源码

def _mvee(self, points, tolerance):
        # Taken from: http://stackoverflow.com/questions/14016898/port-matlab-bounding-ellipsoid-code-to-python
        points = np.asmatrix(points)
        n, d = points.shape
        q = np.column_stack((points, np.ones(n))).T
        err = tolerance + 1.0
        u = np.ones(n) / n
        while err > tolerance:
            # 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))