在有限点集的凸包中找到最小点

日期2009-12-02 (最后修改), 2009-11-30 (创建)

基于 Philip Wolf [1] 的工作和 Kazuyuki Sekitani 和 Yoshitsugu Yamamoto [2] 的递归算法。

[2] 中的算法有 3 个 epsilon 来避免算法三个部分中的比较问题。下面的代码有一些更改，只有一个 epsilon。更改的目的是避免无限循环。

代码¶

在 [1] 中

from numpy import array, matrix, sin, sqrt, dot, cos, ix_, zeros, concatenate, abs, log10, exp, ones
from numpy.linalg import norm

from mpmath import mpf, mp
mp.dps=80

def find_min_point(P):
#    print "Calling find_min with P: ", P

    if len(P) == 1:
        return P[0]

    eps = mpf(10)**-40

    P = [array([mpf(i) for i in p]) for p in P]

    # Step 0. Choose a point from C(P)
    x  = P[array([dot(p,p) for p in P]).argmin()]

    while True:

        # Step 1. \alpha_k := min{x_{k-1}^T p | p \in P}
        p_alpha = P[array([dot(x,p) for p in P]).argmin()]

        if dot(x,x-p_alpha) < eps:
            return array([float(i) for i in x])

        Pk = [p for p in P if abs(dot(x,p-p_alpha)) < eps]

        # Step 2. P_k := { p | p \in P and x_{k-1}^T p = \alpha_k}
        P_Pk = [p for p in P if not array([(p == q).all() for q in Pk]).any()]

        if len(Pk) == len(P):
            return array([float(i) for i in x])

        y = find_min_point(Pk)


        p_beta = P_Pk[array([dot(y,p) for p in P_Pk]).argmin()]

        if dot(y,y-p_beta) < eps:
            return array([float(i) for i in y])


        # Step 4. 
        P_aux = [p for p in P_Pk if (dot(y-x,y-p)>eps) and (dot(x,y-p)!=0)]
        p_lambda = P_aux[array([dot(y,y-p)/dot(x,y-p) for p in P_aux]).argmin()]
        lam = dot(x,p_lambda-y) / dot(y-x,y-p_lambda)

        x += lam * (y-x)



if __name__ == '__main__':
    print find_min_point( [array([ -4.83907292e+00,   2.22438863e+04,  -2.67496763e+04]), array([   9.71147604, -351.46404195, -292.18064276]), array([  4.60452808e+00,   1.07020174e+05,  -1.25310230e+05]), array([  2.16080134e+00,   5.12019937e+04,  -5.96167833e+04]), array([  2.65472146e+00,   6.70546443e+04,  -7.71619656e+04]), array([  1.55775358e+00,  -1.34347516e+05,   1.53209265e+05]), array([   13.22464295,  1869.01251292, -2137.61850989])])


    print find_min_point( [array([ -4.83907292e+00,   2.22438863e+04,  -2.67496763e+04]), array([   9.71147604, -351.46404195, -292.18064276]), array([  4.60452808e+00,   1.07020174e+05,  -1.25310230e+05]), array([  2.16080134e+00,   5.12019937e+04,  -5.96167833e+04]), array([  2.65472146e+00,   6.70546443e+04,  -7.71619656e+04]), array([  1.55775358e+00,  -1.34347516e+05,   1.53209265e+05]), array([   13.22464295,  1869.01251292, -2137.61850989]), array([ 12273.18670123,  -1233.32015854,  61690.10864825])])

[ 13.0643029   -3.03446491  -2.65980139]
[  1.61870596e-04  -3.78774039e-05  -3.29329552e-05]

参考文献¶

1. 在多面体中找到最近的点
2. 一种递归算法，用于在多面体中找到最小范数点，以及两个多面体中的一对最近点

章节作者: JorgeEduardoCardona