分母很小时如何计算数字

Question

我尝试用两个随机点的多面体方法来寻找最优点。 但我发现我得到了不准确的数字，因为我需要将一些数字除以非常小的零。 那我有什么办法处理这个问题吗？ 我试图在分母中加1，但它给了我错误的答案，特别是对于Rosenbrock函数。 这是我的密码。你知道吗

"""
Created on Wed Nov  1 11:50:28 2017

@author: su
"""
import numpy as np
def sorting(f,x,y,z):
    A = np.array([[x[0], x[1], f(x[0], x[1])],
     [y[0], y[1], f(y[0], y[1])],
     [z[0], z[1], f(z[0], z[1])]
     ])
    A = A[A[:,2].argsort()]
    x = A[2][:2]
    y = A[1][:2] 
    z = A[0][:2]
    return x,y,z
def pmin(f,x,y,z,epsilon):
    x,y,z = sorting(f,x,y,z)
    cond1 = float(np.linalg.norm(y-z)) / float(1+np.linalg.norm(z))
    cond2 = float(np.linalg.norm(x-z)) / float(1+np.linalg.norm(z))
    p = np.zeros((1,2))
    count = 0
    while cond1 > epsilon or cond2 > epsilon or count < 10000:
        p = findp(x,y,z)
        if f(x[0],x[1]) > f(p[0],p[1]) :
            x = p
        elif f(y[0],y[1]) > f(p[0],p[1]) :
            y = p
        elif f(z[0],z[1]) > f(p[0],p[1]) :
            z = p
        else:
            x, y = shrink(z,x,y)
            #x, y = shrink(z,y,x)
        x,y,z = sorting(f,x,y,z)
        #print(x,y,z)
        cond1 = float(np.linalg.norm(y-z)) / float(1+np.linalg.norm(z))
        cond2 = float(np.linalg.norm(x-z)) / float(1+np.linalg.norm(z))
        count += 1
    return p, 

def findp(x,y,z):
    p = y + float(np.dot((x-y),(z-y)))/float(np.dot((z-y),(z-y))) * (z-y)
    #p = y + float(np.dot((x-y),(z-y)))/float(1+np.dot((z-y),(z-y))) * (z-y) to avoid float division by zero
    return p-(x-p)
def shrink(x,y,z):
    new_y = x+0.5*(y-x)
    new_z = x+0.5*(z-x)
    return new_y, new_z
def f1(x1,x2):
    return 2*(x1-4)**2 + 3*(x2-5)**2

def f2(x1,x2):
    return (1-x1)**2 + 100*(x2-x1**2)**2

points = np.array([[1,2],[3,4],[9,10]])
eps = float(1e-10)
min1 = pmin(f1,points[0],points[1],points[2],eps)
min2 = pmin(f2,points[0],points[1],points[2],eps)
print "point1"
print "minimum from first function:", min1
print "minimum from second function:", min2
points = np.array([[0,5],[1,3],[1,0]])
min1 = pmin(f1,points[0],points[1],points[2],eps)
min2 = pmin(f2,points[0],points[1],points[2],eps)
print "point2"
print "minimum from first function:", min1
print "minimum from second function:", min2
print(min2)
#from scipy.optimize import fmin
#print fmin(f1, x0=np.array([0]),args=(0,))
#print fmin(f2, x0=np.array([0]),args=(0,))

Answer 1

问题是命名器和分母都很小，甚至精确到0.0，之后算法崩溃。试试这个版本的findp，它应该在数值上更稳定：

def findp(x, y, z):
  nom = float(np.dot((x - y), (z - y)))
  denom = float(np.dot((z - y), (z - y)))
  ratio = 1e-10 if nom < 1e-15 and denom < 1e-15 else nom / denom
  p = y + ratio * (z - y)
  return p - (x - p)