对于多项式，获取其所有极值并通过突出显示所有单调部分来绘制它
For a polynomial, get all its extrema and plot it by highlighting all monotonic pieces
有人问我这个有趣的问题，我认为值得将其张贴在这里，因为 Stack Overflow 上没有任何相关主题。
假设我在长度 n 向量 pc 中有多项式系数，其中变量 x 的度数 n - 1 的多项式可以用其原始形式表示:
pc[1] + pc[2] * x + pc[3] * x ^ 2 + ... + pc[n] * x ^ (n - 1)
R核心函数polyroot可以求出这个多项式在复数域中的所有根。但通常我们也对极值感兴趣，对于单变量函数，局部最小值和最大值交替出现，将函数分解成单调的片段。
我的问题是：
- 如何求多项式实域的所有极值（实际上是所有鞍点）？
- 如何用两种配色方案绘制此多项式：红色代表上升部分，绿色代表下降部分？
最好将其写成一个函数，以便我们可以轻松探索/可视化多项式。
例如，考虑一个 5 次多项式：
pc <- c(1, -2.2, -13.4, -5.1, 1.9, 0.52)
获取多项式的所有鞍点
事实上，鞍点可以通过对多项式的一阶导数使用polyroot来找到。这是一个函数。
SaddlePoly <- function (pc) {
## a polynomial needs be at least quadratic to have saddle points
if (length(pc) < 3L) {
message("A polynomial needs be at least quadratic to have saddle points!")
return(numeric(0))
}
## polynomial coefficient of the 1st derivative
pc1 <- pc[-1] * seq_len(length(pc) - 1)
## roots in complex domain
croots <- polyroot(pc1)
## retain roots in real domain
## be careful when testing 0 for floating point numbers
rroots <- Re(croots)[abs(Im(croots)) < 1e-14]
## note that `polyroot` returns multiple root with multiplicies
## return unique real roots (in ascending order)
sort(unique(rroots))
}
xs <- SaddlePoly(pc)
#[1] -3.77435640 -1.20748286 -0.08654384 2.14530617
计算多项式
我们需要能够计算多项式才能绘制它。 定义了一个函数 g 可以计算多项式及其任意导数。这里我把这个函数复制进去重命名为PolyVal.
PolyVal <- function (x, pc, nderiv = 0L) {
## check missing aruments
if (missing(x) || missing(pc)) stop ("arguments missing with no default!")
## polynomial order p
p <- length(pc) - 1L
## number of derivatives
n <- nderiv
## earlier return?
if (n > p) return(rep.int(0, length(x)))
## polynomial basis from degree 0 to degree `(p - n)`
X <- outer(x, 0:(p - n), FUN = "^")
## initial coefficients
## the additional `+ 1L` is because R vector starts from index 1 not 0
beta <- pc[n:p + 1L]
## factorial multiplier
beta <- beta * factorial(n:p) / factorial(0:(p - n))
## matrix vector multiplication
base::c(X %*% beta)
}
例如，我们可以计算多项式在其所有鞍点处的值：
PolyVal(xs, pc)
#[1] 79.912753 -4.197986 1.093443 -51.871351
为单调件画出具有 2 色方案的多项式
这是一个查看/探索多项式的函数。
ViewPoly <- function (pc, extend = 0.1) {
## get saddle points
xs <- SaddlePoly(pc)
## number of saddle points (if 0 the whole polynomial is monotonic)
n_saddles <- length(xs)
if (n_saddles == 0L) {
message("the polynomial is monotonic; program exits!")
return(NULL)
}
## set a reasonable xlim to include all saddle points
if (n_saddles == 1L) xlim <- c(xs - 1, xs + 1)
else xlim <- extendrange(xs, range(xs), extend)
x <- c(xlim[1], xs, xlim[2])
## number of monotonic pieces
k <- length(xs) + 1L
## monotonicity (positive for ascending and negative for descending)
y <- PolyVal(x, pc)
mono <- diff(y)
ylim <- range(y)
## colour setting (red for ascending and green for descending)
colour <- rep.int(3, k)
colour[mono > 0] <- 2
## loop through pieces and plot the polynomial
plot(x, y, type = "n", xlim = xlim, ylim = ylim)
i <- 1L
while (i <= k) {
## an evaluation grid between x[i] and x[i + 1]
xg <- seq.int(x[i], x[i + 1L], length.out = 20)
yg <- PolyVal(xg, pc)
lines(xg, yg, col = colour[i])
i <- i + 1L
}
## add saddle points
points(xs, y[2:k], pch = 19)
## return (x, y)
list(x = x, y = y)
}
我们可以通过以下方式可视化问题中的示例多项式：
ViewPoly(pc)
#$x
#[1] -4.07033952 -3.77435640 -1.20748286 -0.08654384 2.14530617 2.44128930
#
#$y
#[1] 72.424185 79.912753 -4.197986 1.093443 -51.871351 -45.856876
另一种解决方案，使用 R 包 polynom.
重新实现 SaddlePoly 和 PolyVal
library(polynom)
SaddlePoly <- function (pc) {
## a polynomial needs be at least quadratic to have saddle points
if (length(pc) < 3L) {
message("A polynomial needs be at least quadratic to have saddle points!")
return(numeric(0))
}
## polynomial coefficient of the 1st derivative
## pc1 <- pc[-1] * seq_len(length(pc) - 1) ## <- removed
## roots in complex domain
croots <- solve(deriv(polynomial(pc))) ## <- use package "polynom"
## retain roots in real domain
## be careful when testing 0 for floating point numbers
rroots <- Re(croots)[abs(Im(croots)) < 1e-14]
## note that `polyroot` returns multiple root with multiplicies
## return unique real roots (in ascending order)
sort(unique(rroots))
}
xs <- SaddlePoly(pc)
#[1] -3.77435640 -1.20748286 -0.08654384 2.14530617
## a complete re-implementation using package "polynom"
PolyVal <- function (x, pc, nderiv = 0L) {
## check missing aruments
if (missing(x) || missing(pc)) stop ("arguments missing with no default!")
## create "polynomial" object
p <- polynomial(pc)
## take derivatives
i <- 0
while (i < nderiv) {
p <- deriv(p)
i <- i + 1
}
## evaluate "polynomial" with "predict"
predict(p, x)
}
PolyVal(xs, pc)
#[1] 79.912753 -4.197986 1.093443 -51.871351
## use `ViewPoly` as it is
ViewPoly(pc)
#$x
#[1] -4.07033952 -3.77435640 -1.20748286 -0.08654384 2.14530617 2.44128930
#
#$y
#[1] 72.424185 79.912753 -4.197986 1.093443 -51.871351 -45.856876
在我看来，polynom 包使多项式的构造变得容易。 poly.calc 函数允许从其根或拉格朗日插值构造多项式。
## (x - 1) ^ 3
p1 <- poly.calc(rep(1, 3))
## x * (x - 1) * (x - 2) * (x - 3)
p2 <- poly.calc(0:3)
## Lagrange interpolation through 0:4 and rnorm(5, 0:4, 1)
set.seed(0); x <- 0:4; y <- rnorm(5, 0:4, 1)
p3 <- poly.calc(x, y)
要查看这些多项式，我们可以使用 polynom 或 PolyView 中的函数 plot.polynomial。但是，这两个函数在选择 xlim 作为绘图时具有不同的逻辑。
par(mfrow = c(3, 2), mar = c(4, 4, 1, 1))
## plot `p1`
plot(p1)
ViewPoly(unclass(p1))
## plot `p2`
plot(p2)
ViewPoly(unclass(p2))
## plot `p3`
plot(p3)
ViewPoly(unclass(p3))
有人问我这个有趣的问题，我认为值得将其张贴在这里，因为 Stack Overflow 上没有任何相关主题。
假设我在长度 n 向量 pc 中有多项式系数，其中变量 x 的度数 n - 1 的多项式可以用其原始形式表示:
pc[1] + pc[2] * x + pc[3] * x ^ 2 + ... + pc[n] * x ^ (n - 1)
R核心函数polyroot可以求出这个多项式在复数域中的所有根。但通常我们也对极值感兴趣，对于单变量函数，局部最小值和最大值交替出现，将函数分解成单调的片段。
我的问题是：
- 如何求多项式实域的所有极值（实际上是所有鞍点）？
- 如何用两种配色方案绘制此多项式：红色代表上升部分，绿色代表下降部分？
最好将其写成一个函数，以便我们可以轻松探索/可视化多项式。
例如，考虑一个 5 次多项式：
pc <- c(1, -2.2, -13.4, -5.1, 1.9, 0.52)
获取多项式的所有鞍点
事实上，鞍点可以通过对多项式的一阶导数使用polyroot来找到。这是一个函数。
SaddlePoly <- function (pc) {
## a polynomial needs be at least quadratic to have saddle points
if (length(pc) < 3L) {
message("A polynomial needs be at least quadratic to have saddle points!")
return(numeric(0))
}
## polynomial coefficient of the 1st derivative
pc1 <- pc[-1] * seq_len(length(pc) - 1)
## roots in complex domain
croots <- polyroot(pc1)
## retain roots in real domain
## be careful when testing 0 for floating point numbers
rroots <- Re(croots)[abs(Im(croots)) < 1e-14]
## note that `polyroot` returns multiple root with multiplicies
## return unique real roots (in ascending order)
sort(unique(rroots))
}
xs <- SaddlePoly(pc)
#[1] -3.77435640 -1.20748286 -0.08654384 2.14530617
计算多项式
我们需要能够计算多项式才能绘制它。 g 可以计算多项式及其任意导数。这里我把这个函数复制进去重命名为PolyVal.
PolyVal <- function (x, pc, nderiv = 0L) {
## check missing aruments
if (missing(x) || missing(pc)) stop ("arguments missing with no default!")
## polynomial order p
p <- length(pc) - 1L
## number of derivatives
n <- nderiv
## earlier return?
if (n > p) return(rep.int(0, length(x)))
## polynomial basis from degree 0 to degree `(p - n)`
X <- outer(x, 0:(p - n), FUN = "^")
## initial coefficients
## the additional `+ 1L` is because R vector starts from index 1 not 0
beta <- pc[n:p + 1L]
## factorial multiplier
beta <- beta * factorial(n:p) / factorial(0:(p - n))
## matrix vector multiplication
base::c(X %*% beta)
}
例如，我们可以计算多项式在其所有鞍点处的值：
PolyVal(xs, pc)
#[1] 79.912753 -4.197986 1.093443 -51.871351
为单调件画出具有 2 色方案的多项式
这是一个查看/探索多项式的函数。
ViewPoly <- function (pc, extend = 0.1) {
## get saddle points
xs <- SaddlePoly(pc)
## number of saddle points (if 0 the whole polynomial is monotonic)
n_saddles <- length(xs)
if (n_saddles == 0L) {
message("the polynomial is monotonic; program exits!")
return(NULL)
}
## set a reasonable xlim to include all saddle points
if (n_saddles == 1L) xlim <- c(xs - 1, xs + 1)
else xlim <- extendrange(xs, range(xs), extend)
x <- c(xlim[1], xs, xlim[2])
## number of monotonic pieces
k <- length(xs) + 1L
## monotonicity (positive for ascending and negative for descending)
y <- PolyVal(x, pc)
mono <- diff(y)
ylim <- range(y)
## colour setting (red for ascending and green for descending)
colour <- rep.int(3, k)
colour[mono > 0] <- 2
## loop through pieces and plot the polynomial
plot(x, y, type = "n", xlim = xlim, ylim = ylim)
i <- 1L
while (i <= k) {
## an evaluation grid between x[i] and x[i + 1]
xg <- seq.int(x[i], x[i + 1L], length.out = 20)
yg <- PolyVal(xg, pc)
lines(xg, yg, col = colour[i])
i <- i + 1L
}
## add saddle points
points(xs, y[2:k], pch = 19)
## return (x, y)
list(x = x, y = y)
}
我们可以通过以下方式可视化问题中的示例多项式：
ViewPoly(pc)
#$x
#[1] -4.07033952 -3.77435640 -1.20748286 -0.08654384 2.14530617 2.44128930
#
#$y
#[1] 72.424185 79.912753 -4.197986 1.093443 -51.871351 -45.856876
另一种解决方案，使用 R 包 polynom.
SaddlePoly 和 PolyVal
library(polynom)
SaddlePoly <- function (pc) {
## a polynomial needs be at least quadratic to have saddle points
if (length(pc) < 3L) {
message("A polynomial needs be at least quadratic to have saddle points!")
return(numeric(0))
}
## polynomial coefficient of the 1st derivative
## pc1 <- pc[-1] * seq_len(length(pc) - 1) ## <- removed
## roots in complex domain
croots <- solve(deriv(polynomial(pc))) ## <- use package "polynom"
## retain roots in real domain
## be careful when testing 0 for floating point numbers
rroots <- Re(croots)[abs(Im(croots)) < 1e-14]
## note that `polyroot` returns multiple root with multiplicies
## return unique real roots (in ascending order)
sort(unique(rroots))
}
xs <- SaddlePoly(pc)
#[1] -3.77435640 -1.20748286 -0.08654384 2.14530617
## a complete re-implementation using package "polynom"
PolyVal <- function (x, pc, nderiv = 0L) {
## check missing aruments
if (missing(x) || missing(pc)) stop ("arguments missing with no default!")
## create "polynomial" object
p <- polynomial(pc)
## take derivatives
i <- 0
while (i < nderiv) {
p <- deriv(p)
i <- i + 1
}
## evaluate "polynomial" with "predict"
predict(p, x)
}
PolyVal(xs, pc)
#[1] 79.912753 -4.197986 1.093443 -51.871351
## use `ViewPoly` as it is
ViewPoly(pc)
#$x
#[1] -4.07033952 -3.77435640 -1.20748286 -0.08654384 2.14530617 2.44128930
#
#$y
#[1] 72.424185 79.912753 -4.197986 1.093443 -51.871351 -45.856876
在我看来，polynom 包使多项式的构造变得容易。 poly.calc 函数允许从其根或拉格朗日插值构造多项式。
## (x - 1) ^ 3
p1 <- poly.calc(rep(1, 3))
## x * (x - 1) * (x - 2) * (x - 3)
p2 <- poly.calc(0:3)
## Lagrange interpolation through 0:4 and rnorm(5, 0:4, 1)
set.seed(0); x <- 0:4; y <- rnorm(5, 0:4, 1)
p3 <- poly.calc(x, y)
要查看这些多项式，我们可以使用 polynom 或 PolyView 中的函数 plot.polynomial。但是，这两个函数在选择 xlim 作为绘图时具有不同的逻辑。
par(mfrow = c(3, 2), mar = c(4, 4, 1, 1))
## plot `p1`
plot(p1)
ViewPoly(unclass(p1))
## plot `p2`
plot(p2)
ViewPoly(unclass(p2))
## plot `p3`
plot(p3)
ViewPoly(unclass(p3))