Правильная реализация scipy.optimize для решения всех параметров для плоского уравнения

Question

Трехмерное уравнение плоскости определяется как:

a * x + b * y + c * z + d = 0

где x, y и z - координаты точек данных, а a, b, c и d определяют параметры плоскости. У меня есть несколько точек данных x, y, z, которые, как я знаю, примерно образуют плоскость. Теперь я хочу подогнать эти точки данных к плоскости и найти параметры a, b, c, d. У меня пока есть эта функция:

import scipy.optimize as op
def _plane(ws,x,y,z):
    cost = np.sum(ws[0]*x + ws[1]*y + ws[2]*z + ws[3])
    print cost, np.sum(ws)
    return cost
out = op.minimize(_plane,ws,args=(x,y,z),method='SLSQP',options={'maxiter':1000, 'disp':1})

ws = plane([np.mean(x),np.mean(y),np.mean(z),0.001],x,y,z)

Выход:

130.78467 -0.3011288588643074
130.78467 -0.3011288588643074
130.78467 -0.3011288439631462
130.78467 -0.3011288439631462
130.78468 -0.3011288439631462
130.78467 -0.3011288439631462
-119765.375 -1024.3011288588643
-119765.375 -1024.3011288588643
-119765.375 -1024.3011288439632
-119765.375 -1024.3011288439632
-119765.375 -1024.3011288439632
-119765.375 -1024.3011288439632
Optimization terminated successfully.    (Exit mode 0)
            Current function value: -119765.375
            Iterations: 2
            Function evaluations: 12
            Gradient evaluations: 2
out:      fun: -119765.375
     jac: array([0., 0., 0., 0.])
 message: 'Optimization terminated successfully.'
    nfev: 12
     nit: 2
    njev: 2
  status: 0
 success: True
       x: array([-4.07030255e-01, -8.12448636e-02, -1.02381385e+03,  1.00000000e-03])

Как видите, cost всегда уменьшается, и значение d никогда не меняется, независимо от того, что я даю изначально. Является ли эта реализация правильной?

UPDATE

Использование cost = abs(np.sum(ws[0]*x + ws[1]*y + ws[2]*z + ws[3])) решает проблему оптимизации, а установка method='cg' дает наилучшие результаты.

Answer 1

Для решения конкретной задачи настройки плоскости на набор точек можно использовать более простой подход:

>>> import numpy as np
>>> points = np.array([(3,0,0), (0,4,0), (3,0,1),
...                    (0,4,1), (6,-4,0), (6,-4,1)])

Вычислить центроид и вычесть его из каждой точки.

>>> centroid = np.average(points, axis=0)
>>> centered = points - centroid

Рассчитать разложение по одному значению массива центрированных точек.Третий собственный вектор нормален к плоскости и имеет норму 1, так что вы можете напрямую использовать его в качестве коэффициентов a, b, c плоскости.

>>> _, _, vh = np.linalg.svd(centered)
>>> a, b, c = normal = vh[2]

Центроид должен лежать на плоскости, поэтомупоследний коэффициент (d) может быть вычислен из вектора нормалей и центроида.

>>> d = -normal @ centroid
>>> plane_coefficients = a, b, c, d
>>> print(plane_coefficients)
(0.8, 0.6, 0.0, -2.4000000000000004)