R 使用啥类型的正交多项式？

Posted 2023-03-12

【中文标题】R 使用啥类型的正交多项式？

【英文标题】：What type of orthogonal polynomials does R use?R 使用什么类型的正交多项式？

【发布时间】：2018-05-17 14:37:56

【问题描述】：

我试图匹配 R 中以下代码中的正交多项式：

X <- cbind(1, poly(x = x, degree = 9))

但是在 python 中。

为此，我实现了自己的给出正交多项式的方法：

def get_hermite_poly(x,degree):
    #scipy.special.hermite()
    N, = x.shape
    ##
    X = np.zeros( (N,degree+1) )
    for n in range(N):
        for deg in range(degree+1):
            X[n,deg] = hermite( n=deg, z=float(x[deg]) )
    return X

虽然它似乎不匹配它。有人知道它使用的正交多项式类型吗？我尝试在文档中搜索但没有说。

---

为了提供一些上下文，我正在尝试在 python (https://stats.stackexchange.com/questions/313265/issue-with-convergence-with-sgd-with-function-approximation-using-polynomial-lin/315185#comment602020_315185) 中实现以下 R 代码：

set.seed(1234)

N <- 10
x <- seq(from = 0, to = 1, length = N)
mu <- sin(2 * pi * x * 4)
y <- mu
plot(x,y)

X <- cbind(1, poly(x = x, degree = 9))
# X <- sapply(0:9, function(i) x^i)
w <- rnorm(10)

learning_rate <- function(t) .1 / t^(.6)

n_samp <- 2
for(t in 1:100000) 
  mu_hat <- X %*% w
  idx <- sample(1:N, n_samp)
  X_batch <- X[idx,]
  y_batch <- y[idx]
  score_vec <- t(X_batch) %*% (y_batch - X_batch %*% w)

  change <- score_vec * learning_rate(t)
  w <- w + change


plot(mu_hat, ylim = c(-1, 1))
lines(mu)
fit_exact <- predict(lm(y ~ X - 1))
lines(fit_exact, col = 'red')
abs(w - coef(lm(y ~ X - 1)))

因为它似乎是唯一一种使用梯度下降和多项式特征的线性回归的方法。

我觉得任何正交多项式（或至少正交）都应该工作并给出条件编号为 1 的粗麻布，但我似乎无法让它在 python 中工作。相关问题：How does one use Hermite polynomials with Stochastic Gradient Descent (SGD)?

【参考方案1】：

poly 使用 QR 分解，如 this answer 中的一些详细描述。

我认为您真正想要的是如何使用 python 复制 R 的 poly 的输出。

在这里，我编写了一个基于 R 实现的函数。我还添加了一些 cmets，以便您可以看到 R 中的等效语句是什么样的：

import numpy as np

def poly(x, degree):
    xbar = np.mean(x)
    x = x - xbar

    # R: outer(x, 0L:degree, "^")
    X = x[:, None] ** np.arange(0, degree+1)

    #R: qr(X)$qr
    q, r = np.linalg.qr(X)

    #R: r * (row(r) == col(r))
    z = np.diag((np.diagonal(r)))  

    # R: Z = qr.qy(QR, z)
    Zq, Zr = np.linalg.qr(q)
    Z = np.matmul(Zq, z)

    # R: colSums(Z^2)
    norm1 = (Z**2).sum(0)

    #R: (colSums(x * Z^2)/norm2 + xbar)[1L:degree]
    alpha = ((x[:, None] * (Z**2)).sum(0) / norm1 +xbar)[0:degree]

    # R: c(1, norm2)
    norm2 = np.append(1, norm1)

    # R: Z/rep(sqrt(norm1), each = length(x))
    Z = Z / np.reshape(np.repeat(norm1**(1/2.0), repeats = x.size), (-1, x.size), order='F')

    #R: Z[, -1]
    Z = np.delete(Z, 0, axis=1)
    return [Z, alpha, norm2];

检查这是否有效：

x = np.arange(10) + 1
degree = 9
poly(x, degree)

返回矩阵的第一行是

[-0.49543369,  0.52223297, -0.45342519,  0.33658092, -0.21483446,
  0.11677484, -0.05269379,  0.01869894, -0.00453516],

与R中的相同操作相比

poly(1:10, 9)
# [1] -0.495433694  0.522232968 -0.453425193  0.336580916 -0.214834462
# [6]  0.116774842 -0.052693786  0.018698940 -0.004535159

【讨论】：

哦，哇，你链接的帖子太长了。就像一个简单的问题一样，我认为使用像 qr 这样的因式分解会人为地降低多项式的次数。我记得当我的学位比数据多时进行了分解，它将我的数据矩阵X 减少到数据点的数量。你的方法也是这样吗？我真的不希望它那样做，没关系，我知道我有太多的度数，我是故意这样做的，但即使我有太多的度数/参数，我希望多项式仍然是正交的（不改变我的数字单项式是至关重要的）。

学位不能多于数据。不确定我是否真的了解您的用例如何需要它。但这可能是一个比这里更适合在stats.stackexchange.com 上讨论的问题。

不能是什么意思？你当然可以。只需选择更高次的多项式。我不是要你评估它的统计可靠性。我在问您建议的方法是否会影响我使用的单项式的数量。谢谢你的时间顺便说一句！ :)

我并不是建议您应该使用哪种特定方法。你的问题是“poly 在 R 中使用什么方法”和“我怎样才能在 Python 中得到同样的东西”？我试图展示。 R的算法可以比数据有更多的度数吗？试试poly(1:10, degree=11) 并亲自看看。有关哪种方法最适合您的特定应用程序的相关（但不同）问题的更详细建议，您将在交叉验证时获得更多运气。这个网站更适合编程方面的问题，但对于数学的基本问题，你会得到更多专业人士的关注。

numpy.linalg.qr 和 scipy.linalg.qr 都连接到相同的 LAPACK 例程，所以我不确定为什么会有差异。但是，如果 numpy 函数对您更有效，我会更新答案以改用它。