2

Here is the source code of R poly function (boundary checking are removed). Why we can use QR to build polynomial expansion, which is very different from orthogonal polynomial tutorials like this one. Could anyone help me to understand the connection between these two?

xbar <- mean(x)
x <- x - xbar
X <- outer(x, 0L:degree, "^")
QR <- qr(X)

z <- QR$qr
z <- z * (row(z) == col(z))
Z <- qr.qy(QR, z)
norm2 <- colSums(Z^2)
alpha <- (colSums(x * Z^2)/norm2 + xbar)[1L:degree]
norm2 <- c(1, norm2)

Z <- Z/rep(sqrt(norm2[-1L]), each = length(x))
colnames(Z) <- 0L:degree
Z <- Z[, -1, drop = FALSE]
Haitao Du
  • 32,885
  • 17
  • 118
  • 213
  • 1
    Could you expand your question to be independent of the code presented? QR decomposition is equivalent to Gram Schmidt orthogonalization, in which case it builds a sequence of orthogonal polynomials that approximate your function with minimal least-squares error. – Alex R. Jul 22 '16 at 19:21
  • 2
    I think there is an interesting question here: How is orthogonalization in a function space related to orthogonalization in of the projection of that space into $R^n$ induced by evaluation at a collection of data points. – Matthew Drury Jul 22 '16 at 19:33
  • Thanks @AlexR. I think Matt may be better to help me to expand my question. Basically I read what is QR decomposition from linear algebra book, and read what is orthogonal polynomial, but not quite understand why R can use such code to give me the answer for "orthogonal polynomial expansion". I am not able to explain clearly is because my knowledge is not sufficient to describe this in math and want some help to see the math behind these lines of R code. – Haitao Du Jul 22 '16 at 19:40
  • @hxd1011 Is that what you were after? Do you want me to take a shot at editing your question? – Matthew Drury Jul 22 '16 at 23:47
  • @MatthewDrury we discussed this question couple days ago, if you have time, could you please help me to revise the question? – Haitao Du Aug 05 '16 at 15:27

0 Answers0