Check if point lies inside mulit-dimensional ellipsoid

Question

I have problem with determining if certain points lie inside multi-dimensional ellipsoid. The real problem is caused by potential transformations (rotations, translations) and I don't know how to apply them to my computations.

In my work I will have set of points for which I need to find minimum volume enclosing ellipsoid. Then I will have to check other set of points and determine which of them belong to the found ellipsoid.

I have decided to use code from here:

function [] = Checkup()
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]];
P = points\'; % <- remove the \ symbol here
dimen = length(P(:,1));

% c - vector with ellipse centers
[A, c] = MinVolEllipse(P, 0.01);
[~, Q, V] = svd(A);
radiuses = 1:dimen;

% Calculate radiuses
for i = 1:dimen
    radiuses(1, i) = 1 / sqrt(Q(i,i));
end

% Check if points lie within ellipse and print Ok for every point inside ellipsoid
for i = 1:length(P(1,:)) % length(P(1,:)) is number of points
    value = 0;
    for j = 1:dimen
        %adding ((p_i - c_i) / (r_i))^2 value
        value = value + ( ( (P(j,i) - c(j, 1) )^2 ) / (radiuses(1, j)^2));

    end
    if value <= 1
        disp('Ok')
    end
end

Right now this code does not print Ok text (but it should print it 8 times).

According to this: "V is the rotation matrix that gives you the orientation of the ellipsoid" I think this is the last element I need to use for my code to work.

My approach

Ok so I made some changes to my code. Right now I try to do something like this:

Assuming I have center of ellipsoid and its radius.

For every point:

1. Update position as follows: point = point - center -> in other words translate it as if ellipsoid had it center in point (0, 0, ... 0)

2. Rotate it as if ellipsoid was parallel to all axes: point = invert(V) * point

3. Check if point lies inside ellipsoid using equation: (point.x / radius.x)^2 + ... + (point.i / radius.i)^2 where i is number of dimensions

4. If equation gave result <= 1 then point lies within ellipse.

Is this approach good? I mean - I don't know if I'm allowed to invert V matrix and use it as if it inverted previous rotation...

Proper solution

Looks like my proposed approach was good but there is one even better:

function [result] = Ellipse()
% Returns vector consisting of 10 entries
% that represent error ratio for every test set 

result = 0:9;

% Run tests for all point sets
for pointSet = 0:9
    [Ptraining, Ptest] = LoadPointSet(pointSet);
    count = 0;

    % c - vector with ellipse centers
    [A, c] = MinVolEllipse(Ptraining, 0.00001);

    % Check if points lie within ellipse
    for i = 1:length(Ptest(1,:)) % length(P(1,:)) is number of points
        pp = Ptest(:, i) - c;
        value = pp' * A * pp;    
        if value > 1
            count = count + 1;
        end
    end

    % Get the error ratio
    index = pointSet + 1;
    result(index) = count / length(Ptest(1,:));
end

end

Answer 1

Right now this code does not print Ok text (but it should print it 8 times).

First off, there is no reason to expect that an approximate algorithm such as the MinVolEllipse function will give the exact minimum volume enclosing ellipse. You are supplying a tolerance to MinVolEllipse, so you should expect some error in the result from that function.

More importantly, you don't need to use a singular value decomposition here (but see below for details). The equation for the surface of an ellipsoid is (x-c)^TA(x-c)=1. All you need to check is whether (x-c)^TA(x-c) is less than one (plus some tolerance) for each of your points:

tol_mvee = 0.01;
tol_dist = 0.1;

[A, c] = MinVolEllipse(P, tol_mvee);

% Check if points lie within ellipse and print Ok for every point inside ellipsoid
for i = 1:length(P(1,:)) % length(P(1,:)) is number of points
    x   = P(i,:);
    x_c = x-c;
    d   = dot(x_c, A*x_c);
    if d < 1+tol_dist
        disp('Ok');
    end
end

Your first point is clearly inside the ellipsoid. All of the other seven points are on or very close to the boundary of the true minimum volume enclosing ellipsoid. The algorithm will have a bit of trouble with that, and it will also have a bit of trouble with the fact that the resulting ellipsoid is markedly ellipsoidal rather than spherical (see below on condition number).

A singular value decomposition might be useful to tell you whether the A matrix from MinVolEllipse is ill-conditioned, which is somewhat the case in this particular problem. The condition number for your matrix is over 200, which means you had better not expect to get a result with a tolerance much smaller than 1e-6.