Как найти минимум нелинейной многомерной функции с использованием метода Ньютона (код не линейной алгебры)

Ответ 1

Возможно, вам удастся найти то, что вам нужно, в Numerical Recipes на веб-странице C. Существует свободная версия доступная в Интернете. Здесь (PDF) - это глава, содержащая метод Ньютона-Рафсона, реализованный на C. Вы также можете посмотреть, что доступно в Netlib (LINPack и др.).

Ответ 2

Ссылка на числовые рецепты была наиболее полезной. Я закончил символически дифференцировать мою оценку ошибки, чтобы произвести 30 частных производных, а затем использовал метод Ньютона, чтобы установить их все в ноль. Вот основные моменты кода:

__doc.findzero = [[function(functions, partials, point, [epsilon, steps]) returns table, boolean
Where
  point       is a table mapping variable names to real numbers 
              (a point in N-dimensional space)
  functions   is a list of functions, each of which takes a table like
              point as an argument
  partials    is a list of tables; partials[i].x is the partial derivative
              of functions[i] with respect to 'x'
  epilson     is a number that says how close to zero we're trying to get
  steps       is max number of steps to take (defaults to infinity)
  result      is a table like 'point', boolean that says 'converged'
]]

-- See Numerical Recipes in C, Section 9.6 [http://www.nrbook.com/a/bookcpdf.php]




function findzero(functions, partials, point, epsilon, steps)
  epsilon = epsilon or 1.0e-6
  steps = steps or 1/0
  assert(#functions > 0)
  assert(table.numpairs(partials[1]) == #functions,
         'number of functions not equal to number of variables')
  local equations = { }
  repeat
    if Linf(functions, point) <= epsilon then
      return point, true
    end
    for i = 1, #functions do
      local F = functions[i](point)
      local zero = F
      for x, partial in pairs(partials[i]) do
        zero = zero + lineq.var(x) * partial(point)
      end
      equations[i] = lineq.eqn(zero, 0)
    end
    local delta = table.map(lineq.tonumber, lineq.solve(equations, {}).answers)
    point = table.map(function(v, x) return v + delta[x] end, point)
    steps = steps - 1
  until steps <= 0
  return point, false
end


function Linf(functions, point)
  -- distance using L-infinity norm
  assert(#functions > 0)
  local max = 0
  for i = 1, #functions do
    local z = functions[i](point)
    max = math.max(max, math.abs(z))
  end
  return max
end

Ответ 3

В качестве альтернативы использованию метода Ньютона метод Simplex Nelder-Mead идеально подходит для этой проблемы и упоминается в Numerical Recpies in C.

Rob

Ответ 4

Вы запрашиваете алгоритм минимизации функции. Существует два основных класса: локальный и глобальный. Ваша проблема - наименьшие квадраты, поэтому локальные и глобальные алгоритмы минимизации должны сходиться к одному и тому же уникальному решению. Локальная минимизация намного эффективнее глобальной, поэтому выберите ее.

Существует множество локальных алгоритмов минимизации, но особенно хорошо подходит для задач наименьших квадратов Levenberg-Marquardt. Если у вас нет такого решателя (например, из MINPACK), вы, вероятно, можете уйти с помощью метода Ньютона:

x <- x - (hessian x)^-1 * grad x

где вы вычисляете обратную матрицу, умноженную на вектор, используя линейный решатель.

Ответ 5

Поскольку у вас уже есть частные производные, как насчет общего подхода с градиентом-спуском?

Ответ 6

Возможно, вы думаете, что у вас есть достаточно хорошее решение, но для меня самый простой способ подумать об этом - сначала понять его в случае с 1-переменной, а затем распространить его на матричный случай.

В случае с 1-переменной, если вы разделите первую производную на вторую производную, вы получите (отрицательный) размер шага в следующую пробную точку, например. -V/А.

В случае N-переменной первая производная является вектором, а вторая производная является матрицей (гессиан). Вы умножаете вектор производной на обратную по отношению к второй производной, а результат - отрицательный шаг-вектор на следующую пробную точку, например. -V * (1/А)

Я предполагаю, что вы можете получить матрицу Гессиана 2-й производной. Вам понадобится рутина, чтобы инвертировать ее. Их много в различных пакетах линейной алгебры, и они довольно быстры.

(Для читателей, которые не знакомы с этой идеей, предположим, что две переменные являются 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-я производная), она ответит на один шаг. С другой стороны, если 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);
}

Ответ 7

Мое мнение заключается в использовании стохастического оптимизатора, например метода частиц Рога.