Wie findet man das Minimum einer nichtlinearen, multivariaten Funktion mit der Newton-Methode (Code, nicht lineare Algebra)?

Question

Ich versuche, eine Parameterschätzung durchzuführen und möchte Parameterschätzungen wählen, die den quadratischen Fehler in einer vorhergesagten Gleichung minimieren über etwa 30 Variablen . Wäre die Gleichung linear, würde ich einfach die 30 partiellen Ableitungen berechnen, sie alle auf Null setzen und einen Löser für lineare Gleichungen verwenden. Aber leider die Gleichung ist nichtlinear und ihre Ableitungen sind es auch.

Wäre die Gleichung über eine einzige Variable, würde ich einfach Folgendes verwenden Newton-Verfahren (auch bekannt als Newton-Raphson). Das Internet ist reich an Beispielen und Code zur Umsetzung der Newton-Methode für Funktionen mit einer einzigen Variablen .

Da ich etwa 30 Variablen habe, wie kann ich eine numerische Lösung für dieses Problem mit der Newton-Methode programmieren? ? Ich habe die Gleichung in geschlossener Form und kann die erste und zweite Ableitung berechnen, aber ich weiß nicht so recht, wie es weitergehen soll. Ich habe im Internet eine große Anzahl von Behandlungen gefunden, aber sie gehen schnell in die schwere Matrixnotation über. Ich habe gefunden etwas mäßig Hilfreiches auf Wikipedia, aber ich habe Probleme, es in Code zu übersetzen.

Ich befürchte, dass ich bei der Matrixalgebra und den Matrixinvertierungen nicht weiterkomme. Ich kann eine Matrix mit einem Löser für lineare Gleichungen invertieren, aber ich mache mir Sorgen um die richtigen Zeilen und Spalten, die Vermeidung von Transpositionsfehlern und so weiter.

Um ganz konkret zu sein:

• Ich möchte mit Tabellen arbeiten, die Variablen auf ihre Werte abbilden. Ich kann eine Funktion für eine solche Tabelle schreiben, die den quadratischen Fehler bei einer solchen Tabelle als Argument zurückgibt. Ich kann auch Funktionen erstellen, die eine partielle Ableitung in Bezug auf eine bestimmte Variable zurückgeben.

• Ich habe eine vernünftige Ausgangsschätzung für die Werte in der Tabelle, daher mache ich mir keine Sorgen um die Konvergenz.

• Ich bin mir nicht sicher, wie ich die Schleife schreiben soll, die eine Schätzung (Wertetabelle für jede Variable), die Funktion und eine Tabelle mit partiell abgeleiteten Funktionen verwendet, um eine neue Schätzung zu erstellen.

Bei letzterem bräuchte ich Hilfe. Jede direkte Hilfe oder Hinweise auf gute Quellen werden sehr geschätzt.

---

Edit: Da ich die erste und zweite Ableitung in geschlossener Form habe, möchte ich sie nutzen und langsamer konvergierende Methoden wie Simplex-Suchen vermeiden.

Answer 1

Der Link zu den numerischen Rezepten war sehr hilfreich. Ich habe meine Fehlerschätzung symbolisch differenziert, um 30 partielle Ableitungen zu erzeugen, und dann die Newton-Methode verwendet, um sie alle auf Null zu setzen. Hier sind die Highlights des Codes:

__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

Answer 2

Vielleicht finden Sie, was Sie brauchen, auf der Webseite Numerical Recipes in C. Dort gibt es eine kostenlose Version online verfügbar . Hier (PDF) ist das Kapitel, in dem die Newton-Raphson-Methode in C implementiert ist. Sie können sich auch ansehen, was unter Netlib (LINPack, et. al.).

Answer 3

Alternativ zur Verwendung der Newton-Methode kann die Simplex-Methode von Nelder-Mead ist ideal für dieses Problem geeignet und wird in Numerical Recpies in C erwähnt.

Rob

Answer 4

Sie fragen nach einem Algorithmus zur Funktionsminimierung. Es gibt zwei Hauptklassen: lokal und global. Da es sich bei Ihrem Problem um die kleinsten Quadrate handelt, sollten sowohl lokale als auch globale Minimierungsalgorithmen zu derselben eindeutigen Lösung konvergieren. Die lokale Minimierung ist weitaus effizienter als die globale, also wählen Sie diese.

Es gibt viele lokale Minimierungsalgorithmen, aber einer, der sich besonders gut für Probleme der kleinsten Quadrate eignet, ist Levenberg-Marquardt . Wenn Sie keinen solchen Solver zur Hand haben (z. B. von MINPACK), können Sie wahrscheinlich mit der Newton-Methode auskommen:

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

wo Sie die inverse Matrix multipliziert mit einem Vektor unter Verwendung eines linearen Solvers berechnen.

Answer 5

Da Sie bereits über die partiellen Ableitungen verfügen, wie wäre es mit einem allgemeinen Ansatz des Gradientenabstiegs?