12

我正在尝试进行一些参数估计，并希望选择参数估计，以最小化大约 30 个变量的预测方程中的平方误差。如果方程是线性的，我只需计算 30 个偏导数，将它们全部设置为零，然后使用线性方程求解器。但不幸的是，这个方程是非线性的，它的导数也是如此。

如果方程超过一个变量，我会使用牛顿法（也称为 Newton-Raphson）。网络上有丰富的示例和代码来实现单变量函数的牛顿方法。

鉴于我有大约 30 个变量，我如何使用牛顿法对这个问题的数值解决方案进行编程？我有封闭形式的方程，可以计算一阶和二阶导数，但我不知道如何从那里开始。我在网上找到了大量的处理方法，但它们很快就进入了重矩阵符号。我在 Wikipedia 上发现了一些很有帮助的东西，但我在将其翻译成代码时遇到了麻烦。

我担心分解的地方是矩阵代数和矩阵求逆。我可以用线性方程求解器反转矩阵，但我担心得到正确的行和列，避免转置错误等等。

具体来说：

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

  • 我对表中的值有一个合理的起始估计，所以我不担心收敛。

  • 我不确定如何编写使用估计值（每个变量的值表）、函数和偏导函数表来生成新估计值的循环。

最后一个是我想要帮助的。任何直接帮助或指向良好资源的指针都将受到热烈赞赏。

编辑：由于我有封闭形式的一阶和二阶导数，我想利用它们并避免更慢收敛的方法，如单纯形搜索。

4

7 回答 7

4

数字食谱链接最有帮助。我最终象征性地对我的误差估计进行微分以产生 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
于 2009-01-08T20:42:38.480 回答
3

您可能可以在 C 网页中的数值食谱中找到您需要的内容。网上有 免费版本。 这里(PDF) 是包含用 C 实现的 Newton-Raphson 方法的章节。您可能还想查看Netlib中可用的内容（LINPack 等）。

于 2008-12-24T20:34:44.117 回答
3

作为使用牛顿法的替代方法，Nelder-Mead 的单纯形法非常适合解决这个问题，并在 C 中的 Numerical Recpies 中引用。

抢

于 2008-12-24T20:47:51.517 回答
3

您要求函数最小化算法。有两个主要类：本地和全局。您的问题是最小二乘，因此局部和全局最小化算法都应该收敛到相同的唯一解决方案。局部最小化比全局更有效，所以选择它。

有许多局部最小化算法，但一种特别适合最小二乘问题的算法是Levenberg-Marquardt。如果您手头没有这样的求解器（例如来自MINPACK），那么您可能可以摆脱牛顿的方法：

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

您使用线性求解器计算逆矩阵乘以向量。

于 2010-03-16T06:36:58.570 回答
1

既然你已经有了偏导数，那么一般的梯度下降方法怎么样？

于 2009-01-08T15:55:38.530 回答
1

也许你认为你有一个足够好的解决方案，但对我来说，考虑这个问题的最简单方法是先在 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]

这是对称的。）

如果表面是抛物线（常数二阶导数），它将在 1 步中得到答案。另一方面，如果二阶导数非常不恒定，则可能会遇到振荡。将每一步切成两半（或一小部分）应该使其稳定。

如果 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);
}
于 2009-01-14T15:46:19.410 回答
0

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

于 2009-01-14T18:17:05.377 回答