minimization - sistemas - Cómo encontrar un mínimo de función no lineal, multivariante utilizando el método de Newton(código no álgebra lineal)

sistemas de ecuaciones no lineales octave

Question

Estoy tratando de hacer una estimación de parámetros y quiero elegir estimaciones de parámetros que minimicen el error cuadrado en una ecuación pronosticada sobre unas 30 variables . Si la ecuación fuera lineal, solo calcularía las 30 derivadas parciales, las pondría a cero y usaría un solucionador de ecuaciones lineales. Pero desafortunadamente la ecuación es no lineal y también lo son sus derivados.

Si la ecuación fuera sobre una sola variable, solo usaría el método de Newton (también conocido como Newton-Raphson). La Web es rica en ejemplos y códigos para implementar el método de Newton para funciones de una sola variable .

Dado que tengo alrededor de 30 variables, ¿cómo puedo programar una solución numérica para este problema usando el método de Newton ? Tengo la ecuación en forma cerrada y puedo calcular las derivadas primera y segunda, pero no sé muy bien cómo proceder desde allí. He encontrado una gran cantidad de tratamientos en la web, pero rápidamente entran en notación de matriz pesada. He encontrado algo moderadamente útil en Wikipedia, pero tengo problemas para traducirlo en código.

Lo que me preocupa de romper es en el álgebra matricial y las inversiones de la matriz. Puedo invertir una matriz con un solucionador de ecuaciones lineales, pero me preocupa obtener las filas y columnas correctas, evitar los errores de transposición, y así sucesivamente.

Para ser bastante concreto:

• Quiero trabajar con tablas asignando variables a sus valores. Puedo escribir una función de una tabla que devuelve el error cuadrado dado a dicha tabla como argumento. También puedo crear funciones que devuelven una derivada parcial con respecto a cualquier variable dada.

• Tengo una estimación inicial razonable para los valores en la tabla, por lo que no estoy preocupado por la convergencia.

• No estoy seguro de cómo escribir el ciclo que utiliza una estimación (tabla de valores para cada variable), la función y una tabla de funciones derivadas parciales para producir una nueva estimación.

Eso último es con lo que me gustaría ayudar. Se apreciará calurosamente toda ayuda directa o sugerencias de buenas fuentes.

Editar: Como tengo el primero y el segundo derivados en forma cerrada, me gustaría aprovecharlos y evitar los métodos convergentes más lentos, como las búsquedas simples.

Answer 1

Como una alternativa al uso del método de Newton, el Método Simplex de Nelder-Mead es ideal para este problema y se menciona en Numerical Recpies en C.

Robar

Answer 2

Es posible que pueda encontrar lo que necesita en la página web de Recetas Numéricas en C. Hay una versión gratuita disponible en línea . Aquí (PDF) es el capítulo que contiene el método Newton-Raphson implementado en C. También es posible que desee ver lo que está disponible en Netlib (LINPack, et al.).

Answer 3

Como ya tiene las derivadas parciales, ¿qué le parece un enfoque general de descenso de gradiente?

Answer 4

El enlace Recetas Numéricas fue de gran ayuda. Concluí simbólicamente diferenciando mi estimación de error para producir 30 derivadas parciales, luego usé el método de Newton para ponerlos a cero. Estos son los aspectos más destacados del código:

__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 5

Mi opinión es utilizar un optimizador estocástico, por ejemplo, un método de Enjambre de partículas.

Answer 6

Tal vez piense que tiene una solución lo suficientemente buena, pero para mí, la manera más fácil de pensar sobre esto es primero entenderla en el caso de 1 variable y luego extenderla al caso de la matriz.

En el caso de 1 variable, si divide la primera derivada por la segunda derivada, obtendrá el tamaño de paso (negativo) para su próximo punto de prueba, por ejemplo -V / A.

En el caso de la variable N, la primera derivada es un vector y la segunda derivada es una matriz (el hessiano). Multiplica el vector derivado por el inverso de la segunda derivada, y el resultado es el vector paso negativo para su próximo punto de prueba, por ejemplo -V * (1 / A)

Supongo que puedes obtener la segunda matriz Hessiana derivada. Necesitarás una rutina para invertirlo. Hay muchos de estos en varios paquetes de álgebra lineal, y son bastante rápidos.

(Para lectores que no están familiarizados con esta idea, supongamos que las dos variables son x y y, y la superficie es v (x, y). Entonces la primera derivada es el vector:

 V = [ dv/dx, dv/dy ]

y la segunda derivada es la matriz:

 A = [dV/dx]
    [dV/dy]

o:

 A = [ d(dv/dx)/dx, d(dv/dy)/dx]
    [ d(dv/dx)/dy, d(dv/dy)/dy]

o:

 A = [d^2v/dx^2, d^2v/dydx]
    [d^2v/dxdy, d^2v/dy^2]

que es simétrico)

Si la superficie es parabólica (segunda derivada constante), obtendrá la respuesta en 1 paso. Por otro lado, si la segunda derivada es muy no constante, podrías encontrar oscilación. Cortar cada paso por la mitad (o alguna fracción) debería hacerlo estable.

Si N == 1, verá que hace lo mismo que en el caso de 1 variable.

Buena suerte.

Agregado: Usted quería código:

 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);
}

Answer 7

Usted está solicitando un algoritmo de minimización de funciones. Hay dos clases principales: local y global. Su problema son los mínimos cuadrados, por lo que los algoritmos de minimización tanto locales como globales deben converger a la misma solución única. La minimización local es mucho más eficiente que la global, así que seleccione eso.

Hay muchos algoritmos de minimización locales, pero uno particularmente adecuado para problemas de mínimos cuadrados es Levenberg-Marquardt . Si no tiene a mano un solucionador como este (por ejemplo, de MINPACK), entonces probablemente pueda salirse con la suya con el método de Newton:

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

donde se calcula la matriz inversa multiplicada por un vector utilizando un solucionador lineal.