Численно точное линейное программирование в python для проверки возможности линейного разделения точек? - PullRequest
Новая
       40

Численно точное линейное программирование в python для проверки возможности линейного разделения точек?

0 голосов
/ 06 мая 2018

Я использую линейную программу для проверки, разделяются ли два набора точек линией. Теоретически это работает, но на практике возникают проблемы с числовыми значениями, если точки близки к коллинеарным или если они близки друг к другу.

Я использую пакет 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 требует от меня ввода не менее 3 баллов )

  2. Это ошибка в scipy linprog? Или я просто не правильно его использую?

1 Ответ

0 голосов
/ 07 мая 2018

Вот очищенная версия, которая работает. Я удалил пороговую переменную, так как она может быть поглощена знаком до вектора нормалей. Недостатком является то, что мы должны проверить два случая. Я чувствую, что скачок в нуле в пороговой переменной отбросил решатель. 

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
...