python - python中用于检查点是否可以线性分离的数值精确线性规划?
问题描述
我正在使用线性程序来测试两组点是否可以被一条线分开。理论上这是可行的,但在实践中,如果这些点接近共线或彼此靠近,则似乎存在数值问题。
我使用的包是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
我的问题:
如何可靠(有效)解决检测两组点是否线性可分的问题?(另一种方法是构建凸包,我大部分都写出来了。使用 Qhull 有一些问题:Qhull Convex hull 希望我输入至少 3 个点)
这是 scipy linprog 中的错误吗?还是我没有正确使用它?
解决方案
这是一个似乎可以工作的清理版本。我删除了阈值变量,因为它可以吸收到法线向量中。缺点是我们必须检查两种情况。我的直觉是阈值变量中的零跳跃使求解器失效。
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
推荐阅读
- python - 识别元组中的元素
- c# - Adding records to database with many-to-many relations in ASP.NET Core app with EF Core
- c# - C# XML如何停止&从转换为&
- bash - 我的 systemd 单元文件和 bash 脚本不适用于接口 ppp0 检查
- pytorch-lightning - 在 DDP Pytorch Lightning 中跨 GPU 拆分训练数据
- erlang - Erlang - 如何将原子转换为用'“'包裹?
- macos - 如何将 Homebrew 安装到 /opt/?
- amazon-web-services - AWS Glue:NotADirectoryError: [Errno 20] 不是目录:'/tmp/glue-python-scripts-XXXXX/.YYYY' -> '/tmp/glue-python-scripts-XXXX/'
- vue.js - 如何限制按钮的增量?
- python - 分配可以放入计算机的张量的资源耗尽错误