首页 文章

如何使用牛顿方法(代码非线性代数)找到最小的非线性,多元函数

提问于
浏览
12

我正在尝试进行一些参数估计,并希望选择最小化预测方程 over about 30 variables 中的平方误差的参数估计 . 如果方程是线性的,我只计算30个偏导数,将它们全部设为零,并使用线性方程求解器 . 但不幸的是 the equation is nonlinear 以及它的衍生品 .

如果方程式超过单个变量,我只会使用Newton's method(也称为Newton-Raphson) . Web上有丰富的示例和代码来实现Newton的单个变量函数的方法 .

鉴于我有大约30个变量, how can I program a numeric solution to this problem using Newton's method ?我有封闭形式的等式,可以计算一阶和二阶导数,但我没有在维基百科上找到something moderately helpful,但我在将其转换为代码时遇到了麻烦 .

我担心崩溃的地方是矩阵代数和矩阵求逆 . 我可以使用线性方程求解器反转矩阵,但我担心得到正确的行和列,避免换位错误,等等 .

要非常具体:

  • 我想使用表将变量映射到它们的值 . 我可以写一个这样一个表的函数,它返回给出这样一个表作为参数的平方误差 . 我还可以创建相对于任何给定变量返回偏导数的函数 .

  • 我对表中的值有一个合理的初始估计,所以我不担心收敛 .

  • 我不确定如何编写使用估计的循环(每个变量的值表),函数和部分导数函数表来产生新的估计 .

最后一点是我要帮助的 . 任何直接的帮助或指向好的来源将受到热烈的赞赏 .


编辑:因为我有封闭形式的第一和第二衍生物,我想利用它们并避免更简单的融合方法,如单面搜索 .

7 回答

  • 0

    既然你已经有了偏导数,那么一般的梯度下降方法呢?

  • 3

    您可以在C网页的Numerical Recipes中找到所需内容 . 有free version available online . Here(PDF)是包含用C实现的Newton-Raphson方法的章节 . 您可能还想查看Netlib(LINPack等人)提供的内容 .

  • 3

    Numerical Recipes链接最有帮助 . 我结束了象征性地区分我的误差估计以产生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
    
  • 1

    作为使用牛顿方法的替代方法,Simplex Method of Nelder-Mead非常适合这个问题,并在C中的数值重合中引用 .

  • 1

    您要求的是函数最小化算法 . 有两个主要类:本地和全球 . 您的问题是最小二乘,因此局部和全局最小化算法应该收敛到相同的唯一解 . 局部最小化比全局最有效,所以选择它 .

    有许多局部最小化算法,但一个特别适合最小二乘问题的算法是Levenberg-Marquardt . 如果你没有't have such a solver to hand (e.g. from MINPACK) then you can probably get away with Newton'方法:

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

    使用线性求解器计算逆矩阵乘以向量的位置 .

  • 3

    也许你认为你有一个足够好的解决方案,但对我来说,最简单的思考方法是首先在1变量的情况下理解它,然后将其扩展到矩阵的情况 .

    在1变量的情况下,如果将一阶导数除以二阶导数,则得到(负)步长到下一个试验点,例如: -V / A .

    在N变量情况下,一阶导数是矢量,二阶导数是矩阵(Hessian) . 您将导数向量乘以二阶导数的倒数,结果是到下一个试验点的负步进向量,例如, -V *(1 / A)

    我假设你可以得到二阶导数的Hessian矩阵 . 你需要一个例程来反转它 . 在各种线性代数包中有很多这些,它们非常快 .

    (对于不熟悉这个想法的读者,假设两个变量是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]
    

    这是对称的 . )

    如果表面是抛物线的(常数二阶导数),它将在一步得到答案 . 另一方面,如果二阶导数非常不恒定,则可能会遇到振荡 . 将每一步切成两半(或一些部分)应使其稳定 .

    如果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);
    }
    
  • 4

    我的意见是使用随机优化器,例如粒子群方法 .

相关问题