Может быть, вы думаете, что у вас есть достаточно хорошее решение, но для меня самый простой способ подумать об этом - сначала понять его в случае с 1 переменной, а затем распространить его на случай матрицы.
В случае с 1 переменной, если вы разделите первую производную на вторую производную, вы получите (отрицательный) размер шага до следующей пробной точки, например, -V / A.
В случае N-переменной первая производная является вектором, а вторая производная является матрицей (гессианом). Вы умножаете вектор производной на значение, обратное второй производной, и в результате получается отрицательный шаговый вектор к следующей точке испытания, например, -V * (1 / А)
Полагаю, вы можете получить матрицу Гессиана со второй производной. Вам понадобится рутина, чтобы инвертировать его. Их много в различных пакетах линейной алгебры, и они довольно быстрые.
(Для читателей, которые не знакомы с этой идеей, предположим, что две переменные - это x и y, а поверхность - v (x, y). Тогда первой производной является вектор:
V = [ dv/dx, dv/dy ]
и вторая производная является матрицей:
A = [dV/dx]
[dV/dy]
или
A = [ d(dv/dx)/dx, d(dv/dy)/dx]
[ d(dv/dx)/dy, d(dv/dy)/dy]
или
A = [d^2v/dx^2, d^2v/dydx]
[d^2v/dxdy, d^2v/dy^2]
что симметрично.)
Если поверхность параболическая (постоянная 2-я производная), она получит ответ за 1 шаг. С другой стороны, если 2-я производная очень не постоянна, вы можете столкнуться с колебаниями. Разрезание каждого шага пополам (или некоторой доли) должно сделать его стабильным.
Если N == 1, вы увидите, что он делает то же самое, что и в случае с 1 переменной.
Удачи.
Добавлено: Вы хотели код:
double X[N];
// Set X to initial estimate
while(!done){
double V[N]; // 1st derivative "velocity" vector
double A[N*N]; // 2nd derivative "acceleration" matrix
double A1[N*N]; // inverse of A
double S[N]; // step vector
CalculateFirstDerivative(V, X);
CalculateSecondDerivative(A, X);
// A1 = 1/A
GetMatrixInverse(A, A1);
// S = V*(1/A)
VectorTimesMatrix(V, A1, S);
// if S is small enough, stop
// X -= S
VectorMinusVector(X, S, X);
}