The problem can be solved by finding a feasible 8 point of a Linear Program. If you're interested 7 in the full details, as opposed to just 6 plugging an LP into an existing solver, I'd 5 recommend reading Chapter 11.4 in Boyd and Vandenberghe's excellent book on convex optimization.

Set A = (X[1] X[2] ... X[n]), that 4 is, the first column is v1, the second v2, etc.

Solve 3 the following LP problem,

minimize (over x): 1
s.t.     Ax = P
         x^T * [1] = 1
         x[i] >= 0  \forall i

where

1. x^T is the transpose of x
2. [1] is the all-1 vector.

The problem 2 has a solution iff the point is in the convex 1 hull.

The point lies outside of the convex hull 12 of the other points if and only if the direction 11 of all the vectors from it to those other 10 points are on less than one half of a circle/sphere/hypersphere 9 around it.

Here is a sketch for the situation 8 of two points, a blue one inside the convex 7 hull (green) and the red one outside:

For 6 the red one, there exist bisections of the 5 circle, such that the vectors from the point 4 to the points on the convex hull intersect 3 only one half of the circle. For the blue 2 point, it is not possible to find such a 1 bisection.

You don't have to compute convex hull itself, as 11 it seems quite troublesome in multidimensional 10 spaces. There's a well-known property of convex hulls:

Any vector (point) v inside 9 convex hull of points [v1, v2, .., vn] can be presented 8 as sum(ki*vi), where 0 <= ki <= 1 and sum(ki) = 1. Correspondingly, no point 7 outside of convex hull will have such representation.

In 6 m-dimensional space, this will give us the 5 set of m linear equations with n unknowns.

edit
I'm 4 not sure about complexity of this new problem 3 in general case, but for m = 2 it seems linear. Perhaps, somebody 2 with more experience in this area will correct 1 me.

I had the same problem with 16 dimensions. Since 21 even qhull didn't work properly as too much 20 faces had to be generated, I developed my 19 own approach by testing, whether a separating 18 hyperplane can be found between the new 17 point and the reference data (I call this 16 "HyperHull" ;) ).

The problem 15 of finding a separating hyperplane can be 14 transformed to a convex quadratic programming 13 problem (see: SVM). I did this in python using 12 cvxopt with less then 170 lines of code 11 (including I/O). The algorithm works without 10 modification in any dimension even if there 9 exists the problem, that as higher the dimension 8 as higher the number of points on the hull 7 (see: On the convex hull of random points in a polytope). Since the hull isn't explicitely 6 constructed but only checked, whether a 5 point is inside or not, the algorithm has 4 very big advantages in higher dimensions 3 compared to e.g. quick hull.

This algorithm 2 can 'naturally' be parallelized and speed 1 up should be equal to number of processors.

Though the original post was three years 27 ago, perhaps this answer will still be of 26 help. The Gilbert-Johnson-Keerthi (GJK) algorithm 25 finds the shortest distance between two 24 convex polytopes, each of which is defined 23 as the convex hull of a set of generators---notably, the 22 convex hull itself does not have to be calculated. In 21 a special case, which is the case being 20 asked about, one of the polytopes is just 19 a point. Why not try using the GJK algorithm 18 to calculate the distance between P and 17 the convex hull of the points X? If that 16 distance is 0, then P is inside X (or at 15 least on its boundary). A GJK implementation 14 in Octave/Matlab, called ClosestPointInConvexPolytopeGJK.m, along 13 with supporting code, is available at http://www.99main.com/~centore/MunsellAndKubelkaMunkToolbox/MunsellAndKubelkaMunkToolbox.html . A 12 simple description of the GJK algorithm 11 is available in Sect. 2 of a paper, at http://www.99main.com/~centore/ColourSciencePapers/GJKinConstrainedLeastSquares.pdf . I've 10 used the GJK algorithm for some very small 9 sets X in 31-dimensional space, and had 8 good results. How the performance of GJK 7 compares to the linear programming methods 6 that others are recommending is uncertain 5 (although any comparisons would be interesting). The 4 GJK method does avoid computing the convex 3 hull, or expressing the hull in terms of 2 linear inequalities, both of which might 1 be time-consuming. Hope this answer helps.

Are you willing to accept a heuristic answer 14 that should usually work but is not guaranteed 13 to? If you are then you could try this 12 random idea.

Let f(x) be the cube of the 11 distance to P times the number of things 10 in X, minus the sum of the cubes of the 9 distance to all of the points in X. Start 8 somewhere random, and use a hill climbing 7 algorithm to maximize f(x) for x in a sphere 6 that is very far away from P. Excepting 5 degenerate cases, if P is not in the convex 4 hull this should have a very good probability 3 of finding the normal to a hyperplane which 2 P is on one side of, and everything in X 1 is on the other side of.

A write-up to test if a point is in a hull space, using scipy.optimize.minimize.

Based on user1071136's answer.

It does go 6 a lot faster if you compute the convex hull, so 5 I added a couple of lines for people who 4 want to do that. I switched from graham 3 scan (2D only) to the scipy qhull algorithm.

scipy.optimize.minimize 2 documentation:
https://docs.scipy.org/doc/scipy/reference/optimize.nonlin.html

import numpy as np
import scipy.optimize
import matplotlib.pyplot as plt
from scipy.spatial import ConvexHull


def hull_test(P, X, use_hull=True, verbose=True, hull_tolerance=1e-5, return_hull=True):
    if use_hull:
        hull = ConvexHull(X)
        X = X[hull.vertices]

    n_points = len(X)

    def F(x, X, P):
        return np.linalg.norm( np.dot( x.T, X ) - P )

    bnds = [[0, None]]*n_points # coefficients for each point must be > 0
    cons = ( {'type': 'eq', 'fun': lambda x: np.sum(x)-1} ) # Sum of coefficients must equal 1
    x0 = np.ones((n_points,1))/n_points # starting coefficients

    result = scipy.optimize.minimize(F, x0, args=(X, P), bounds=bnds, constraints=cons)

    if result.fun < hull_tolerance:
        hull_result = True
    else:
        hull_result = False

    if verbose:
        print( '# boundary points:', n_points)
        print( 'x.T * X - P:', F(result.x,X,P) )
        if hull_result: 
            print( 'Point P is in the hull space of X')
        else: 
            print( 'Point P is NOT in the hull space of X')

    if return_hull:
        return hull_result, X
    else:
        return hull_result

Test on some sample data:

n_dim = 3
n_points = 20
np.random.seed(0)

P = np.random.random(size=(1,n_dim))
X = np.random.random(size=(n_points,n_dim))

_, X_hull = hull_test(P, X, use_hull=True, hull_tolerance=1e-5, return_hull=True)

Output:

# boundary points: 14
x.T * X - P: 2.13984259782e-06
Point P is in the hull space of X

Visualize 1 it:

rows = max(1,n_dim-1)
cols = rows
plt.figure(figsize=(rows*3,cols*3))
for row in range(rows):
    for col in range(row, cols):
        col += 1
        plt.subplot(cols,rows,row*rows+col)
        plt.scatter(P[:,row],P[:,col],label='P',s=300)
        plt.scatter(X[:,row],X[:,col],label='X',alpha=0.5)
        plt.scatter(X_hull[:,row],X_hull[:,col],label='X_hull')
        plt.xlabel('x{}'.format(row))
        plt.ylabel('x{}'.format(col))
plt.tight_layout()