Newton-Raphson Method

求解非线性方程/方程组

求解非线性方程组，即寻找能够使 \(f(\boldsymbol{X})=\boldsymbol{0}\) 的 \(\boldsymbol{X}\).

具体实现上，有多种算法，如不动点迭代法（fixed-point iteration）、牛顿法（Newton-Raphson method）、割线法（secant method）、二分法（bisection method）等。不动点法的收敛速度较慢，也更可能导致不收敛；最常用的是牛顿法和二分法。

手写 Newton-Raphson 算法

一元方程

Newton_root <- function(f, start,
                        tolerance = 1e-7, max_iter = 50) {
  root <- start
  for (i in 1:max_iter) {
    f_prime <- numDeriv::genD(f, root)$D[1, 1]
    # $D 是一个矩阵，其中不仅包含一阶导数，还有二阶导数
    dx <- f(root) / f_prime
    root <- root - dx
    print(str_c("root = ", root))
    if (abs(dx) < tolerance) {
      return(list(
        root = root, f.root = f(root),
        iterator = i, estimate_precise = abs(dx)
      ))
    }
  }
  print("Maximum number of iterations exceeded.")
}

transformation0 <- function(x) {
  x^3 - sin(x)^2
}

Newton_root(f = transformation0, start = 1)

#> [1] "root = 0.860369147925258"
#> [1] "root = 0.809891223391771"
#> [1] "root = 0.802931615098361"
#> [1] "root = 0.802803774492958"
#> [1] "root = 0.802803731737894"

#> $root
#> [1] 0.8028037
#> 
#> $f.root
#> [1] 4.440892e-15
#> 
#> $iterator
#> [1] 5
#> 
#> $estimate_precise
#> [1] 4.275506e-08

多元方程组

Newton_Raphson <- function(transformation, start, tolerance = 1e-7, max_iter = 50) {
  dimension <- length(start)
  root <- start
  for (i in 1:max_iter) {
    f <- transformation(root)
    Jacobi <- numDeriv::genD(transformation, root)$D[1:dimension, 1:dimension]
    dx <- solve(Jacobi, f)
    root <- root - dx
    str_c(root, collapse = ", ") %>%
      str_c("root = [", ., "]") %>%
      print()

    precise <- sqrt(sum(dx * dx))
    if (precise < tolerance) {
      return(list(
        root = root, f.root = transformation(root),
        iterator = i, estimate_precise = precise
      ))
    }
  }
  print("Maximum number of iterations exceeded.")
}

transformation1 <- function(x) {
  c(
    f1 = x[1] + x[2] - 2.5,
    f2 = x[1] * x[2] - 1
  )
}

Newton_Raphson(transformation1, start = c(0, 3))

#> [1] "root = [0.33333333333384, 2.16666666666507]"
#> [1] "root = [0.484848484848594, 2.01515151514307]"
#> [1] "root = [0.499849984998516, 2.00015001500144]"
#> [1] "root = [0.499999984999999, 2.000000015]"
#> [1] "root = [0.5, 2]"

#> $root
#> [1] 0.5 2.0
#> 
#> $f.root
#>            f1            f2 
#>  0.000000e+00 -2.220446e-16 
#> 
#> $iterator
#> [1] 5
#> 
#> $estimate_precise
#> [1] 2.121321e-08

rootSolve 包实现 Newton-Raphson 算法

算法详情见 CRAN 上的 rootSolve.pdf，可能对标准的 Newton-Raphson 算法进行了若干优化，防止发散

一元方程

寻找 root，使传入的函数作用于 root 时返回的标量为0

fs <- function(s) {
  s^3 - 4 * s^2 - 10 * s + 4
}

# 一次性在4个初始位置开始迭代，求出了所有三个数值解
rootSolve::multiroot(f = fs, start = c(-5, 0, 2.5, 8))

#> $root
#> [1] -2.0000000  0.3542487  0.3542487  5.6457513
#> 
#> $f.root
#> [1] -6.375274e-10 -8.881784e-16  4.440892e-16  1.421085e-14
#> 
#> $iter
#> [1] 7
#> 
#> $estim.precis
#> [1] 1.593857e-10

# 求出参数区间内的所有解
rootSolve::uniroot.all(f = fs, interval = c(-10, 10))

#> [1] -2.000000  0.354230  5.645754

二元方程组 without Jacobi

一般来说，不需要显示给出偏导数的 Jacobi 矩阵

transformation2 <- function(x) {
  c(
    f1 = 10 * x[1] + 3 * x[2]^2 - 3,
    f2 = x[1]^2 - exp(x[2]) - 2
  )
}
rootSolve::multiroot(f = transformation2, start = c(1, 1))

#> $root
#> [1] -1.445552 -2.412158
#> 
#> $f.root
#>            f1            f2 
#>  5.117684e-12 -6.084022e-14 
#> 
#> $iter
#> [1] 10
#> 
#> $estim.precis
#> [1] 2.589262e-12

二元方程组 with Jacobi

显示给出偏导数的 Jacobi 矩阵，在某些时候可以加快迭代速度

derivs <- function(x) {
  c(
    10, 6 * x[2], # f11, f12
    2 * x[1], -exp(x[2]) # f21, f22
  ) %>%
    matrix(nrow = 2, byrow = T)
}
rootSolve::multiroot(f = transformation2, start = c(0, 0), jacfunc = derivs)

#> $root
#> [1] -1.445552 -2.412158
#> 
#> $f.root
#>            f1            f2 
#>  1.166651e-09 -1.390243e-11 
#> 
#> $iter
#> [1] 29
#> 
#> $estim.precis
#> [1] 5.902766e-10

方程组中有若干参数

transformation3 <- function(x, parms) {
  c(
    f1 = x[1] + x[2] + x[3]^2 - parms[1],
    f2 = x[1]^2 - x[2] + x[3] - parms[2],
    f3 = 2 * x[1] - x[2]^2 + x[3] - 1
  )
}
parms <- c(12, 2)
rootSolve::multiroot(transformation3, start = c(1, 1, 1), parms = parms)

#> $root
#> [1] 1 2 3
#> 
#> $f.root
#>            f1            f2            f3 
#>  3.087877e-10  4.794444e-09 -8.678146e-09 
#> 
#> $iter
#> [1] 6
#> 
#> $estim.precis
#> [1] 4.593792e-09

rootSolve::multiroot(transformation3, c(0, 0, 0), parms = parms * 2)

#> $root
#> [1] -0.8103629  1.4861978  4.8295098
#> 
#> $f.root
#>            f1            f2            f3 
#> -4.973799e-14  4.143352e-12 -5.280221e-13 
#> 
#> $iter
#> [1] 11
#> 
#> $estim.precis
#> [1] 1.573704e-12

矩阵形式的方程组

rootSolve::multiroot()的第一个参数（非线性变换）可以返回一个矩阵。计算出的解将保证这个矩阵为全零矩阵。

例：求解 25 元非线性方程组

\[ \mathbf{X} \cdot \mathbf{X} \cdot \mathbf{X}=\left[\begin{array}{lllll} 1 & 2 & 3 & 4 & 5 \\ 6 & 7 & 8 & 9 & 10 \\ 11 & 12 & 13 & 14 & 15 \\ 16 & 17 & 18 & 19 & 20 \\ 21 & 22 & 23 & 24 & 25 \end{array}\right] \]

transformation4 <- function(x) {
  X <- matrix(x, nrow = 5)
  X %*% X %*% X - matrix(1:25, nrow = 5, byrow = TRUE)
}

x <- rootSolve::multiroot(transformation4, start = 1:25)$root
x

#>  [1] -0.67506260 -0.03483809  0.60668343  1.24815964  1.88640623 -0.34547775
#>  [7]  0.16864884  0.68350802  1.19974369  1.71457060 -0.01800918  0.37530821
#> [13]  0.76025826  1.15042418  1.54265669  0.31230570  0.57358465  0.84421247
#> [19]  1.10046033  1.36971978  0.64043429  0.77971070  0.91597598  1.06007417
#> [25]  1.19373256

X <- matrix(nrow = 5, x)
X %*% X %*% X

#>      [,1] [,2] [,3] [,4] [,5]
#> [1,]    1    2    3    4    5
#> [2,]    6    7    8    9   10
#> [3,]   11   12   13   14   15
#> [4,]   16   17   18   19   20
#> [5,]   21   22   23   24   25