I had to search around a bit to understand exactly what kind of orthogonality R’s poly function refers to. It seems to be orthogonality with respect to the following inner product:
\langle p, q \rangle = \sum_{i = 1}^n p(x_i) q(x_i) \, ,
where the x_i are the points provided by the user. Let’s label this family of polynomials p_k(x), where k is the degree of the polynomial.
Since you’re only interested in the values p_k(x_i), rather than the coefficients defining the polynomials, what we’re actually trying to calculate is just a matrix Q defined as
Q = \begin{bmatrix}
p_0(x_1) & p_1(x_1) & \cdots & p_m(x_1) \\
\vdots & \vdots & ~ & \vdots \\
p_0(x_n) & p_1(x_n) & \cdots & p_m(x_n)
\end{bmatrix}
which should satisfy the following properties:
- The columns should be orthogonal (this is the orthogonality of the polynomials)
- The column space of Q should be the vector space spanned by polynomials of degree up to m evaluated on the x_i, which is mostly easily expressed as the column space of the Vandermonde matrix V:
V = \begin{bmatrix}
1 & x_1 & x_1^2 & \cdots & x_1^m \\
\vdots & \vdots & \vdots & ~ & \vdots \\
1 & x_n & x_n^2 & \cdots & x_n^m
\end{bmatrix}
- The first k columns of Q should have the same span as the first k columns of V (this ensures that p_k(x) is a polynomial of degree k)
These are exactly the properties satisfied by the factor Q in the “thin” QR decomposition of the Vandermonde matrix, V = Q R. Thus, what we need to do is a) form the Vandermonde matrix, and b) compute its QR decomposition. Each column of the resulting Q contains the values of one of the polynomials p_k(x). Here’s an implementation
using LinearAlgebra
function poly(x, degree)
V = reduce(hcat, (x.^k for k in 0:degree))
F = qr(V)
return F.Q[:, 2:(degree + 1)]
end
Note one last detail: R’s poly doesn’t return the first column, which is just constant; to mirror this behavior, the implementation above drops the first column of Q and returns columns 2:(degree + 1).
Let’s verify that it works:
julia> poly(1:10, 3)
10×3 Matrix{Float64}:
-0.495434 0.522233 0.453425
-0.385337 0.174078 -0.151142
-0.275241 -0.0870388 -0.377854
-0.165145 -0.261116 -0.334671
-0.0550482 -0.348155 -0.12955
0.0550482 -0.348155 0.12955
0.165145 -0.261116 0.334671
0.275241 -0.0870388 0.377854
0.385337 0.174078 0.151142
0.495434 0.522233 -0.453425
The first two columns here match this StackOverflow answer demonstrating usage of R’s poly: https://stackoverflow.com/a/19484326, so I’d say this passes the sniff test.
