# Searching extreme points of polyhedron

In my Uni, my scientific professor asked me to make some researches about the extreme points of polyhedrals. And I did them. I found that there is still no code in public for searching extreme points for polyhedral with n dimensions (n - x's), but polyhedrons are everywhere (CV, game theories, etc.). I wrote a function for this task and made a python library (also there are matrix and array combination maker functions).

All I want is to make this code more optimal and compatible for all python versions (I have some troubles that some times happen when I install it by "pip install lin", but other times no). I want to make the life of people easier and make it more comfortable.

I am asking for you to test this function on your computer and write if you have any bugs, fails or thoughts on how to make it better. I am open to constructive criticism and non-constructive too (it will help to understand if somebody needs it or that is just a waste of time).

All the examples, instructions and code on my GitHub: https://github.com/r4ndompuff/polyhedral_set

import numpy as np
import itertools as it
import math
import re

def permutation(m,n):
return math.factorial(n)/(math.factorial(n-m)*math.factorial(m))

def matrix_combinations(matr,n):
timed = list(map(list, it.combinations(matr, n)))
for i in range(n):
timed[i][i][i] = np.asscalar(timed[i][i][i])
all = np.array(list(timed))
return all

def array_combinations(arr,n):
timed = list(map(list, it.combinations(arr, n)))
for i in range(n):
timed[i][i] = np.asscalar(timed[i][i])
all = np.array(list(timed))
return all

def check_extreme(matr, arr, x, sym_comb, m):
sym_comb = sym_comb.replace(']', '')
sym_comb = sym_comb.replace('[', '')
sym_comb = re.split("[ ,]", sym_comb)
for i in range(m):
if sym_comb[i] == '>':
return 0
elif sym_comb[i] == '>=':
return 0
elif sym_comb[i] == '<':
return 0
elif sym_comb[i] == '<=':
return 0
elif sym_comb[i] == '=':
return 0
elif sym_comb[i] == '!=':
return 0
else:
return 0
return 1

def extreme_points(m,n,A,b,sym_comb):
# Input
A = np.array(A).reshape(m,n)
b = np.array(b).reshape(m,1)
# Proccess
ans_comb = np.zeros((1,n))
arr_comb = array_combinations(b,n)
matr_comb = matrix_combinations(A,n)
for i in range(int(permutation(n,m))):
if np.linalg.det(matr_comb[i]) != 0:
x = np.linalg.solve(matr_comb[i],arr_comb[i])
ans_comb = np.vstack([ans_comb,x])
ans_comb = np.delete(ans_comb, (0), axis=0)
j = 0
for i in range(len(ans_comb)):
if check_extreme(A, b, ans_comb[j], sym_comb, m):
ans_comb = ans_comb
j = j + 1
else:
ans_comb = np.delete(ans_comb, (j), axis=0)
# Output
return ans_comb


• Can you provide some more information on what exactly is going on, from a math standpoint? Is this actually a polyhedron, or a polytope? In how many dimensions? By "extreme", what do you mean? Euclidean norm from (the origin, an arbitrary point)? Commented Apr 21, 2019 at 23:21
• "there is still no code in public for searching extreme points for polyhedral with n dimensions" - I can nearly guarantee that that isn't the case. Commented Apr 21, 2019 at 23:22
• @Reinderien I am still in the process of deciphering the question. I have a math background, and I share the mother tongue with OP. My impression is that by extremal points OP means the vertices of a simplex where a certain linear form (defined by A, b, and the condition) achieves an extremum. I could be wrong.
– vnp
Commented Apr 22, 2019 at 0:32
• @vnp That's kind of what I guessed, and if that's the case, linear programming is a quite well-established field already - with some stuff built right into scipy. Commented Apr 22, 2019 at 0:33
• @Reinderien Agreed. Still deciphering.
– vnp
Commented Apr 22, 2019 at 0:34

I created a rudimentary pull request to your GitHub repo. I won't show all of the content here except for the main file:

import numpy as np
import itertools as it
from math import factorial
import re

def permutation(m, n):
return factorial(n) / (factorial(n - m) * factorial(m))

def matrix_combinations(matr, n):
timed = list(map(list, it.combinations(matr, n)))
return np.array(list(timed))

def array_combinations(arr, n):
timed = list(map(list, it.combinations(arr, n)))
return np.array(list(timed))

def check_extreme(matr, arr, x, sym_comb, m):
sym_comb = sym_comb.replace(']', '')
sym_comb = sym_comb.replace('[', '')
sym_comb = re.split("[ ,]", sym_comb)
for i in range(int(m)):
if sym_comb[i] == '>':
return 0
elif sym_comb[i] == '>=':
return 0
elif sym_comb[i] == '<':
return 0
elif sym_comb[i] == '<=':
return 0
elif sym_comb[i] == '=':
return 0
elif sym_comb[i] == '!=':
return 0
else:
return 0
return 1

def extreme_points(A, b, sym_comb):
# Input
A = np.array(A)
b = np.array(b)
m, n = A.shape
# Process
ans_comb = np.zeros((1, n))
arr_comb = array_combinations(b, n)
matr_comb = matrix_combinations(A, n)
for i in range(int(permutation(n, m))):
if np.linalg.det(matr_comb[i]) != 0:
x = np.linalg.solve(np.array(matr_comb[i], dtype='float'),
np.array(arr_comb[i], dtype='float'))
ans_comb = np.vstack([ans_comb, x])
ans_comb = np.delete(ans_comb, 0, axis=0)
j = 0
for i in range(len(ans_comb)):
if check_extreme(A, b, ans_comb[j], sym_comb, m):
ans_comb = ans_comb
j += 1
else:
ans_comb = np.delete(ans_comb, j, axis=0)
# Output
return ans_comb


Notable changes:

• Do a direct import of factorial
• Don't call asscalar, since it's both unneeded and deprecated
• Don't call a variable all, since that shadows a Python built-in
• Don't need to explicitly pass array dimensions, nor do you need to reshape the arrays
• Drop redundant parens around some expressions
• Use += where applicable
• Fix up almost all PEP8 issues, except for your capital letter A, which is fine in context

This doesn't solve the bigger issue that you should replace 99% of this with a call to scipy. I'll do that separately (I suspect that @vnp is, as well).

• Thank you! In meantime I made it work with any NumPy version now (about loop error you report). Commented Apr 22, 2019 at 7:57
• Sorry, but I can't find the same function in scipy. Can you provide me the name of function? Commented Apr 22, 2019 at 8:17
• @AndreyLovyagin docs.scipy.org/doc/scipy/reference/generated/… Commented Apr 22, 2019 at 16:18
• Thank you very much for all that you have done! Also, I checked scipy.optimize.linprog() function and I can't really understand how it can replace much of the code. I will learn and read more about it, but on first sight, it looks like it can't solve non-equational system (I mean not '=', but any sign, e.g. '>='), that is why I can't see how I can find extreme points by it. Commented Apr 22, 2019 at 20:39
• Read about the A_ub, b_ub arguments. Commented Apr 22, 2019 at 20:40