python中用于检查点是否可以线性分离的数值精确线性规划？

Question

我正在使用线性程序来测试两组点是否可以被一条线分开。理论上这是可行的，但在实践中，如果这些点接近共线或彼此靠近，则似乎存在数值问题。

我使用的包是scipy linprog。

我写了一些代码来举例说明我在说什么。这段代码产生一个包含 N 个点的云“位置”，并选择一些随机线，然后使用这些线周围的边距将“位置”集的子集划分为两个区域 X 和 Y，然后检查线性程序是否成功找到 X 和 Y 之间的分隔线，或者如果线性程序得出结论不存在这样的分隔线（未能找到可行点）。

从输出（下图）可以看出，随着边距从 1 变为 2^(-10)，线性程序正确检测到两个区域是线性可分的概率会衰减。这表明在检测到可以分离非常近的点云时可能存在一些问题。

请注意，因为我正在输入保证线性可分的线性程序集，所以输出应该都是 1。

import numpy as np
from scipy.optimize import linprog

N = 100

for minv in range(10):
    margin = 1/(2**minv)
    locations = np.random.uniform(-1,1,[N,2])
    tests = []
    for i in range(50):
        p = np.random.normal(0,1,[3])
        X = []
        Y = []
        for a in locations:
            if np.dot(a, [p[0],p[1]]) < p[2] - margin:
                X.append(a)
            if np.dot(a, [p[0],p[1]]) > p[2] + margin:
                Y.append(a)
        #X and Y are two linearly seperable subsets of 'locations'
        A = []
        b = []
        if X != [] and Y != []:
            #This if is just to avoid an edge case that causes problems with linprog
            for s in X:
                A.append( [-1*v for v in s] + [1] )
                b.append(-1)
            for s in Y:
                A.append( list(s) + [-1])
                b.append(-1)
            #See: https://www.joyofdata.de/blog/testing-linear-separability-linear-programming-r-glpk/
            res = linprog(c = [0,0,0], A_ub = A, b_ub = b, bounds = [-np.inf, np.inf])
            tests.append(res["success"])
    print(np.mean(tests))

输出：

1.0
1.0
0.909090909091
0.8
0.805555555556
0.5
0.375
0.444444444444
0.378378378378
0.410256410256

我的问题：

1. 如何可靠（有效）解决检测两组点是否线性可分的问题？ （另一种方法是构建凸包，我大部分都写出来了。使用 Qhull 有一些问题：Qhull Convex hull wants me to input at least 3 points）

2. 这是 scipy linprog 中的错误吗？还是我没有正确使用它？

Answer 1

这是一个似乎可以工作的清理版本。我删除了阈值变量，因为它可以吸收到法线向量中。缺点是我们必须检查两种情况。我的直觉是，阈值变量中的零跳跃让求解器崩溃。

import numpy as np
from scipy.optimize import linprog

N = 100

def test(X, Y):
    if not (X.size and Y.size):
        return True, np.full((2,), np.nan)
    A = np.r_[-X, Y]
    b = np.repeat((-1, 1), (X.shape[0], Y.shape[0]))
    res1 = linprog(c=[0,0], A_ub=A, b_ub=b-1e-6, bounds=2*[(None, None)])
    res2 = linprog(c=[0,0], A_ub=A, b_ub=-b-1e-6, bounds=2*[(None, None)])
    if res1['success']:
        return True, res1['x']
    elif res2['success']:
        return True, res2['x']
    else:
        return False, np.full((2,), np.nan)

def split(locations, p, margin):
    proj = np.einsum('ij,j', locations, p[:2])
    X = locations[proj < p[2] - margin]
    Y = locations[proj > p[2] + margin]
    return X, Y

def validate(locations, p, p_recon, margin):
    proj = np.einsum('ij,j', locations, p[:2])
    X = locations[proj < p[2] - margin]
    Y = locations[proj > p[2] + margin]
    if not (X.size and Y.size):
        return True
    return ((np.all(np.einsum('ij,j', X, p_recon) > 1)
             and np.all(np.einsum('ij,j', Y, p_recon) < 1)) 
            or
            (np.all(np.einsum('ij,j', X, p_recon) > -1)
             and np.all(np.einsum('ij,j', Y, p_recon) < -1)))

def random_split(locations):
    inX = np.random.randint(0, 2, (locations.shape[0],), dtype=bool)
    return locations[inX], locations[~inX]

for minv in range(10):
    margin = 1/(2**minv)
    locations = np.random.uniform(-1,1,[N,2])
    P = np.random.normal(0, 1, (50, 3))
    tests, P_recon = zip(*(test(*split(locations, p, margin)) for p in P))
    assert all(validate(locations, *ps, margin) for ps in zip(P, P_recon))
    control = [test(*random_split(locations))[0] for p in P]
    print('test:', np.mean(tests), 'control:', np.mean(control))

样本输出：

test: 1.0 control: 0.0
test: 1.0 control: 0.0
test: 1.0 control: 0.0
test: 1.0 control: 0.0
test: 1.0 control: 0.0
test: 1.0 control: 0.0
test: 1.0 control: 0.0
test: 1.0 control: 0.0
test: 1.0 control: 0.0
test: 1.0 control: 0.0