python 中的数值精确线性规划,用于检查点是否可以线性分离?
Numerically precise linear programming in python for checking whether points can be linearly separated?
我正在使用线性程序来测试两组点是否可以被一条线分开。理论上这是可行的,但在实践中,如果点接近共线或彼此靠近,似乎存在数值问题。
我使用的包是scipy linprog
。
我写了一些代码来举例说明我在说什么。这段代码产生了一个由 N 个点组成的云 'locations',并随机选取了一些线,然后使用这些线周围的边距将 'locations' 集的一个子集划分为两个区域 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 wants me to input at least 3 points)
这是 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
我正在使用线性程序来测试两组点是否可以被一条线分开。理论上这是可行的,但在实践中,如果点接近共线或彼此靠近,似乎存在数值问题。
我使用的包是scipy linprog
。
我写了一些代码来举例说明我在说什么。这段代码产生了一个由 N 个点组成的云 'locations',并随机选取了一些线,然后使用这些线周围的边距将 'locations' 集的一个子集划分为两个区域 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 wants me to input at least 3 points)
这是 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