Linear Polynomial Regression

This tutorial show how to use the P4ML package to perform a naive polynomial regression. This is not really the intended use-case of the package, but it serves to illustrate some basic usage and functionality.

using Polynomials4ML, LinearAlgebra, StaticArrays
using SpheriCart: SphericalHarmonics

Example 1: Univariate Polynomial Regression

First, we specify a target function, which we will try to approximate.

f1(x) = 1 / (1 + 10 * x^2);

Suppose we are given data in the form of argument, function value pairs, with arguments sampled uniformly from the interval [-1, 1].

N = 1_000
Xtrain = 2 * rand(N) .- 1
Ytrain = f1.(Xtrain);

We want to solve the least squares problem

\[ \min_p \sum_{i = 1}^{N} | y_i - p(x_i) |^2\]

for a polynomial $p$. Because $f_1$ is analytic, we use a degree approximately $\sqrt{N}$. Since the distribution of the samples is uniform, the Legrendre polynomials are likely a good choice.

basis = legendre_basis(ceil(Int, sqrt(N)))
c = basis(Xtrain) \ Ytrain
p = x -> sum(c .* basis(x));

A quick statistical test error

Xtest = 2 * rand(N) .- 1
Ytest = f1.(Xtest)
println("Test RMSE: ", norm(Ytest - p.(Xtest)) / sqrt(N))

Test RMSE: 2.685737732226062e-5

Example 2: Polynomial Regression on the Sphere

As a second example we consider a function defined on the unit sphere,

f2(x) = 1 / (1 + 10 * norm(x - SA[1.0,0.0,0.0])^2);

where x is now an SVector{3, Float64} with norm(x) == 1. We generate again uniform samples on the sphere.

N = 1_000
X = [ (x = randn(SVector{3, Float64}); x/norm(x)) for i = 1:N ]
Y = f2.(X);

In this case, spherical harmonics are natural basis functions. These are provided by SpheriCart.jl in both real and complex variants. Since the target function is real, we choose the real basis.

basis = SphericalHarmonics(9)
@show length(basis);

length(basis) = 100

We see that we now have far more basis functions per degree and therefore cannot use as high a degree as before (unless we give ourselves more data). Otherwise we can proceed exactly as above.

c = basis(X) \ Y
p = x -> sum(c .* basis(x));

We can test the RMSE again.

Xtest = [ (x = randn(SVector{3, Float64}); x/norm(x)) for i = 1:N ]
Ytest = f2.(Xtest)
println("Test RMSE: ", norm(Ytest - p.(Xtest)) / sqrt(N));

Test RMSE: 0.009282656369071895

This page was generated using Literate.jl.