Regressione polinomiale con Accord.NET

Questo post mostra un uso della classe PolynomialRegression del framework Accord.NET con l'obiettivo di dimostrare che la regressione polinomiale classica del machine learning può raggiungere interessanti livelli di accuratezza con tempi di learning estremamente brevi.
Anche se le reti neurali di tipo MLP (Multi Layer Perceptron) possono essere considerate degli approssimatori universali di funzioni (si veda Approssimazione con percettroni multistrato altamente configurabili su questo sito web), per alcune tipologie di dataset, un classico algoritmo supervisionato di regressione di machine learning, come quello di regressione polinomiale di cui si tratta qui, può arrivare a livelli di accuratezza accettabili e con un costo computazionale decisamente più basso rispetto a quello di un MLP.

Nel mondo reale i dataset preesistono alla fase di apprendimento, infatti questi vengono ottenuti estraendo dati da database di produzione o file Excel, dall'output di strumenti di misura, dai data-logger collegati a sensori elettronici e così via, e quindi adoperati per le fasi di learning successive; ma poiché qui il focus è l'approssimazione polinomiale in sé e non l'approssimazione di un fenomeno reale, i dataset utilizzati in questo post saranno generati sinteticamente: questo ha il vantaggio di poter stressare l'algoritmo e vedere per quali tipologie di dataset l'algoritmo ha accuratezza accettabile e per quali invece l'algoritmo fatica.

Implementazione

Il codice descritto da questo post è scritto il C#, richiede Microsoft .NET Core 3.1 e il framework Accord.NET 3.8.0.
Per ottenere il codice si veda il paragrafo Download del codice completo in fondo a questo post.
La classe più importante del progetto è la classe PolynomialLeastSquares1to1Regressor che è un wrapper della classe PolynomialRegression di Accord.NET.

using System.Globalization;
using System.Collections.Generic;
using System.Linq;
using Accord.Statistics.Models.Regression.Linear;
using Accord.Math.Optimization.Losses;
using Regressors.Entities;
using Regressors.Exceptions;

namespace Regressors.Wrappers
{
    public class PolynomialLeastSquares1to1Regressor
    {
        private int _degree;
        private bool _isRobust;
        private PolynomialRegression _polynomialRegression;

        public PolynomialLeastSquares1to1Regressor(int degree, bool isRobust = false)
        {
            _degree = degree;
            _isRobust = isRobust;
        }

        public double[] Weights
        {
            get
            {
                AssertAlreadyLearned();
                return _polynomialRegression.Weights;
            }
        }

        public double Intercept
        {
            get
            {
                AssertAlreadyLearned();
                return _polynomialRegression.Intercept;
            }
        }
        private void AssertAlreadyLearned()
        {
            if (_polynomialRegression == null)
                throw new NotTrainedException();
        }

        public string StringfyLearnedPolynomial(string format = "e")
        {
            if (_polynomialRegression == null)
                return string.Empty;
            
            return _polynomialRegression.ToString(format, CultureInfo.InvariantCulture);
        }

        public void Learn(IList<XtoY> dsLearn)
        {            
            double [] inputs = dsLearn.Select(i => i.X).ToArray();
            double [] outputs = dsLearn.Select(i => i.Y).ToArray();
            var pls = new PolynomialLeastSquares() { Degree = _degree, IsRobust = _isRobust };
            _polynomialRegression = pls.Learn(inputs, outputs);
        }

        public IEnumerable<XtoY> Predict(IEnumerable<double> xvalues)
        {
            AssertAlreadyLearned();

            double []xvaluesArray = xvalues.ToArray();
            double [] yvaluesArray = _polynomialRegression.Transform(xvaluesArray);
            for(int i = 0; i < xvaluesArray.Length; ++i)
            {
                yield return new XtoY() {X = xvaluesArray[i], Y = yvaluesArray[i]};
            }
        }

        public static double ComputeError(IEnumerable<XtoY> ds, IEnumerable<XtoY> predicted)
        {
            double [] outputs = ds.Select(i => i.Y).ToArray();
            double [] preds = predicted.Select(i => i.Y).ToArray();
            double error = new SquareLoss(outputs).Loss(preds);
            return error;
        }
    }
}

Esecuzione del codice

Su un sistema Linux, Mac o Windows ove .NET Core 3.1 e git client siano già installati, per scaricare il sorgente e mandare in esecuzione il programma, eseguire i seguenti comandi:

git clone https://github.com/ettoremessina/accord-net-experiments
cd accord-net-experiments/Prototypes/CSharp/Regressors
dotnet build
dotnet run

dotnet run

Started test #PLSR_01
Generating learning dataset
Generating test dataset
Training
Learned polynomial y(x) = 5.000000e-001x^3 + -2.000000e+000x^2 + -3.000000e+000x^1 + -1.000000e+000
Predicting
Error: 9.42921358148529E-24
Saving
Terminated test #PLSR_01

Caso di studio dell'approssimazione polinomiale di una sinusoide

Dato il dataset sintetico generato con una sinusoide $f(x)=\sin x$ nell'intervallo $[-2 \pi, 2 \pi]$, il primo tentativo di approssimazione è effettuato tentando con un polinomio di terzo grado; da questo grafico:

Ragionando empiricamente (attenzione: questo ragionamento non è appricabile in generale), osservando che la funzione da approssimare ha 2 minimi e 2 massimi relativi nell'intervallo considerato ovverosia 5 tratti a direzione di crescita opposta, si osserva che serve almeno un polinomio di grado 5 per realizzare questi 5 tratti a direzione di crescita opposta.
In effetti, tentando ad approssimare con un polinomio di quinto grado le cose iniziano ad andare meglio:

Alzando il grado del polinomio a 10, si nota che ci si avvicina in maniera interessante a una approssimazione molto accurata:

Caso di studio dell'approssimazione polinomiale di una funzione con una cuspide

Sia dato il dataset sintetico generato con la seguente funzione $f(x)=\sqrt{|x|}$ nell'intervallo $[-5, 5]$, che presenta una cuspide nell'origine; la regressione polinomiale ha difficoltà ad approssimare la cuspide come si evince dal seguente grafico che mostra un polinomio di decimo grado che tenta di approssimare la funzione obiettivo:

Anche provando con un polinomio di grado 20:

Caso di studio dell'approssimazione di un tratto di esponenziale

Anche se la funzione esponenziale non è approssimabile con un polinomio, contrariamente a quanto si possa pensare, su un tratto limitato l'approssimazione raggiunge un grande livello di accuratezza. Data la funzione esponenziale $f(x)=e^x$ nell'intervallo $[-5, 5]$, la regressione polinomiale con un polinomio di quinto grado è già accettabile, con un errore quadratico medio pari a 0.9961909800481012, come mostra il grafico:

Utilizzando un polinomio di decimo grado, l'approssimazione diventa molto accurata, con un errore quadratico medio pari a 1.6342296507040001E-06:

Un caso di studio difficile: la sinusoide smorzata

La funzione $f(x)=\frac{\sin 2x}{e^\frac{x}{5}}$, che è una sorta di sinusoide smorzata, nell'intervallo $[-20, 20]$ non sembra facile da approssimare per via dei tanti minimi e massimi locali presenti; per quanto detto sopra si tenta subito ad approssimare con un polinomio di alto grado, ad esempio un ventesimo grado; dal seguente grafico:

Utilizzando un polinomio di quarantesimo grado, l'approssimazione diventa molto accurata, con un errore quadratico medio decisamente più basso pari a 0.03318409290980062:

Alternative

Per esempi di regressione di una funzione reale di una variable con SVM in Weka si veda Regressione con SMO per SVM con kernel PUK in Weka.
Per approfondire invece l'approssimazione di una funzione reale di una variable con un MLP si veda Approssimazione di una funzione reale di una variabile con TensorFlow e Approssimazione di una funzione reale di una variabile con PyTorch.

Download del codice completo

Il codice completo è disponibile su GitHub.
Questo materiale è distribuito su licenza MIT; sentiti libero di usare, condividere, "forkare" e adattare tale materiale come credi.
Sentiti anche libero di pubblicare pull-request e bug-report su questo repository di GitHub oppure di contattarmi sui miei canali social disponibili nell'angolo in alto a destra di questa pagina.