Python中用于三维曲面最小二乘拟合的方法

3

我希望将我的表面方程拟合到一些数据上。我已经尝试了scipy.optimize.leastsq,但是由于我无法指定边界,所以它给了我一个无法使用的结果。我还尝试了scipy.optimize.least_squares,但它给了我一个错误:

ValueError: too many values to unpack

我的方程式是:

f(x,y,z)=(x-A+y-B)/2+sqrt(((x-A-y+B)/2)^2+C*z^2)

需要找到参数 A、B、C,使得在使用以下 x、y、z 点时,上述方程尽可能接近于零:

    [
   [-0.071, -0.85, 0.401],
   [-0.138, -1.111, 0.494],
   [-0.317, -0.317, -0.317],
   [-0.351, -2.048, 0.848]
   ]

这些边界条件为A > 0,B > 0,C > 1。

我应该如何得到这样的拟合?在Python中,什么是最好的工具来完成这项工作?我搜索了关于拟合3D表面的示例,但大多数关于函数拟合的示例都涉及线性或平面拟合。

1个回答

7
我将编辑此答案,提供一个更一般的示例,说明如何使用scipy的通用optimize.minimize方法以及scipy的optimize.least_squares方法来解决此问题。
首先让我们设置问题:
import numpy as np
import scipy.optimize

# ===============================================
# SETUP: define common compoments of the problem


def our_function(coeff, data):
    """
    The function we care to optimize.

    Args:
        coeff (np.ndarray): are the parameters that we care to optimize.
        data (np.ndarray): the input data
    """
    A, B, C = coeff
    x, y, z = data.T
    return (x - A + y - B) / 2 + np.sqrt(((x - A - y + B) / 2) ** 2 + C * z ** 2)


# Define some training data
data = np.array([
    [-0.071, -0.85, 0.401],
    [-0.138, -1.111, 0.494],
    [-0.317, -0.317, -0.317],
    [-0.351, -2.048, 0.848]
])
# Define training target
# This is what we want the target function to be equal to
target = 0

# Make an initial guess as to the parameters
# either a constant or random guess is typically fine
num_coeff = 3
coeff_0 = np.ones(num_coeff)
# coeff_0 = np.random.rand(num_coeff)

这并不是严格的最小二乘法,但可以尝试类似于以下方法:

这种解决方案就像用大锤子敲击问题。可能有一种方法可以使用最小二乘法通过SVD求解器更有效地获得答案,但如果您只是寻找答案,scipy.optimize.minimize将为您找到一个。

# ===============================================
# FORMULATION #1: a general minimization problem

# Here the bounds and error are all specified within the general objective function
def general_objective(coeff, data, target):
    """
    General function that simply returns a value to be minimized.
    The coeff will be modified to minimize whatever the output of this function
    may be.
    """
    # Constraints to keep coeff above 0
    if np.any(coeff < 0):
        # If any constraint is violated return infinity
        return np.inf
    # The function we care about
    prediction = our_function(coeff, data)
    # (optional) L2 regularization to keep coeff small
    # (optional) reg_amount = 0.0
    # (optional) reg = reg_amount * np.sqrt((coeff ** 2).sum())
    losses = (prediction - target) ** 2
    # (optional) losses += reg
    # Return the average squared error
    loss = losses.sum()
    return loss


general_result = scipy.optimize.minimize(general_objective, coeff_0,
                                         method='Nelder-Mead',
                                         args=(data, target))
# Test what the squared error of the returned result is
coeff = general_result.x
general_output = our_function(coeff, data)
print('====================')
print('general_result =\n%s' % (general_result,))
print('---------------------')
print('general_output = %r' % (general_output,))
print('====================')

输出结果如下:
====================
general_result =
 final_simplex: (array([[  2.45700466e-01,   7.93719271e-09,   1.71257109e+00],
       [  2.45692680e-01,   3.31991619e-08,   1.71255150e+00],
       [  2.45726858e-01,   6.52636219e-08,   1.71263360e+00],
       [  2.45713989e-01,   8.06971686e-08,   1.71260234e+00]]), array([ 0.00012404,  0.00012404,  0.00012404,  0.00012404]))
           fun: 0.00012404137498459109
       message: 'Optimization terminated successfully.'
          nfev: 431
           nit: 240
        status: 0
       success: True
             x: array([  2.45700466e-01,   7.93719271e-09,   1.71257109e+00])
---------------------
general_output = array([ 0.00527974, -0.00561568, -0.00719941,  0.00357748])
====================

我在文档中发现,要将其适应于实际的最小二乘法,只需指定计算残差的函数。

# ===============================================
# FORMULATION #2: a special least squares problem

# Here all that is needeed is a function that computes the vector of residuals
# the optimization function takes care of the rest
def least_squares_residuals(coeff, data, target):
    """
    Function that returns the vector of residuals between the predicted values
    and the target value. Here we want each predicted value to be close to zero
    """
    A, B, C = coeff
    x, y, z = data.T
    prediction = our_function(coeff, data)
    vector_of_residuals = (prediction - target)
    return vector_of_residuals


# Here the bounds are specified in the optimization call
bound_gt = np.full(shape=num_coeff, fill_value=0, dtype=np.float)
bound_lt = np.full(shape=num_coeff, fill_value=np.inf, dtype=np.float)
bounds = (bound_gt, bound_lt)

lst_sqrs_result = scipy.optimize.least_squares(least_squares_residuals, coeff_0,
                                               args=(data, target), bounds=bounds)
# Test what the squared error of the returned result is
coeff = lst_sqrs_result.x
lst_sqrs_output = our_function(coeff, data)
print('====================')
print('lst_sqrs_result =\n%s' % (lst_sqrs_result,))
print('---------------------')
print('lst_sqrs_output = %r' % (lst_sqrs_output,))
print('====================')

这里的输出是:
====================
lst_sqrs_result =
 active_mask: array([ 0, -1,  0])
        cost: 6.197329866927735e-05
         fun: array([ 0.00518416, -0.00564099, -0.00710112,  0.00385024])
        grad: array([ -4.61826888e-09,   3.70771396e-03,   1.26659198e-09])
         jac: array([[-0.72611025, -0.27388975,  0.13653112],
       [-0.74479565, -0.25520435,  0.1644325 ],
       [-0.35777232, -0.64222767,  0.11601263],
       [-0.77338046, -0.22661953,  0.27104366]])
     message: '`gtol` termination condition is satisfied.'
        nfev: 13
        njev: 13
  optimality: 4.6182688779976278e-09
      status: 1
     success: True
           x: array([  2.46392438e-01,   5.39025298e-17,   1.71555150e+00])
---------------------
lst_sqrs_output = array([ 0.00518416, -0.00564099, -0.00710112,  0.00385024])
====================

你写的代码给了我正确的结果!不过,我真的不太明白你是如何施加边界条件的? - UN4
哎呀,我忘了做那件事。¯\(ツ)/¯ 我觉得它只是按我的初步猜测使用正数值运行了,并且L2正则化鼓励系数缩小,所以当它达到一个好的答案时就停了下来。然而这种方法非常通用,所以你可以加上 if np.any(coeff < 0): return np.inf 这样就可强制限制条件了。 - Erotemic
啊,好的,那很有道理。谢谢! - UN4
我弄清楚了最小二乘法的实现方式。您指定的函数需要计算残差。在这种情况下,边界可以明确指定。 - Erotemic

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接