 |
 |
 |
Onlangs schreef ik een post over het programmeren van chemometrische algoritmen in het pre-Python-tijdperk, waarbij ik koos voor vintage Fortran als programmeertaal. In dat bericht beloofde ik een algoritme voor Partial Least Squares regressie in een taal uit het pre-Python-tijdperk. In plaats van Fortran presenteer ik dat echter in Algol 68. De voorganger daarvan, Algol 60, werd voornamelijk gebruikt voor onderzoek door computerwetenschappers in de Verenigde Staten en Europa. Als het op serieuze rekenkracht aankwam, konden talen uit de Algol-familie om diverse praktische redenen niet concurreren met Fortran.
In mijn vorige blogpost over dit onderwerp hebben we een eenvoudig OLS-algoritme (Ordinary Least Squares) geprogrammeerd om de lineaire vergelijking A β = y op te lossen voor β, ons “model” in het jargon van machine learning, wanneer er geen exacte oplossing bestaat. In de chemometrie betekenen meetfouten doorgaans dat er geen exacte oplossing is. Als A een dataset zou zijn in de vorm van een m×n-matrix, die m monsters vertegenwoordigt met n waarnemingen per monster, en y de vereiste m uitkomsten voor elk monster zou zijn, dan zou β een kalibratiefactor zijn.
Het berekenen van β is de "training" van ons machine learning "model". In een recent essay heb ik het effect besproken van de antropomorfe taal die in kunstmatige intelligentie wordt gebruikt. Anderen zullen chemometrie associëren met machine learning, maar persoonlijk associeer ik dit prachtige vakgebied met multivariate statistiek en lineaire algebra. Maar dat is mijn persoonlijke mening.
Aangezien datasets in de scheikunde doorgaans gecorreleerde gegevens bevatten, is het beter om de regressie te baseren op de dominante eigenvectoren van de covariantiematrix van je dataset. Lange eigenvectoren duiden op richtingen met grote variatie in een dataset, terwijl korte eigenvectoren wijzen op richtingen met grote correlatie. Vervolgens minimaliseer je de voorspellingsfout naarmate het aantal in aanmerking genomen eigenvectoren toeneemt. Dit proces wordt om voor de hand liggende redenen dimensiereductie genoemd. Een bekende methode is Principal Component Regression (PCR).
In chemische datasets is het waarschijnlijk dat materiaaleigenschappen die voor voorspellingen worden gebruikt, op moleculair niveau gecorreleerd zijn met de te voorspellen eigenschap. Als je bijvoorbeeld spectroscopische metingen gebruikt om materiaaleigenschappen zoals viscositeit te voorspellen, weet je dat beide worden beïnvloed door het gedrag van dezelfde moleculen. Dit is een kanttekening bij PCR, dat de mogelijkheid van een dergelijke correlatie negeert en daardoor dominante eigenvectoren mogelijk verkeerd identificeert. Een betere benadering is Partial Least Squares (PLS)-regressie, die ernaar streeft rekening te houden met die correlatie. PLS vond in de jaren tachtig zijn weg naar de chemometrie en is om redenen die u nu begrijpt een werkpaard geworden.
Laten we eens kijken naar het PLS-algoritme. PLS kan worden getraind om meerdere eigenschappen te voorspellen op basis van één dataset. Vaak wordt echter aangeraden om als eerste stap per gewenste eigenschap een model te trainen en te beoordelen of dat voor jou werkt. Dit wordt het PLS1-algoritme genoemd. Uiteraard kan deze PLS1-code worden uitgevoerd met Algol 68 Genie. Eerst lezen we uit een bestand een dataset met de bekende waarden.
PRAGMAT need gsl PRAGMAT
# Read a training set with (single-column) constituent values.
Training set entries are in a CSV file formatted per record as:
sample name, constituent value, feature values #
INT samples = read int, features = read int;
VOID (open (standin, "training-set.csv", standin channel));
make term (standin, ",");
[samples, features] REAL training set, [samples, 1] REAL constituents;
FOR i TO samples
DO read ((LOC STRING, space, constituents[i, 1], space));
FOR j TO features
DO read (training set[i, j]);
IF j < features
THEN read (space)
FI
OD;
IF i < samples
THEN read (newline)
FI
OD;
Zodra we de gegevens hebben gelezen, definiëren we de vereiste gegevenstypen en centreren we de gegevens rond het gemiddelde.
# Centering columns is required for PCR and PLS, not for OLS. #
MODE MATRIX = FLEX [1 : 0, 1 : 0] REAL,
VECTOR = FLEX [1 : 0] REAL,
COLVEC = MATRIX; # Column vector #
MATRIX col mean set = MEAN training set, COLVEC col mean consts = MEAN constituents;
MATRIX mean set = training set - col mean set,
COLVEC mean consts = constituents - col mean consts;
Hieronder presenteren we een naïeve implementatie van PLS1. We beginnen met het berekenen van de eigenvectoren, gerangschikt in aflopende volgorde van lengte.
COMMENT Partial Least Squares Regression .
NIPALS PLS1 following Ulf Indahl.
COMMENT
MATRIX e := mean set, new set, new consts, eigens,
COLVEC f := mean consts;
TO samples # For this demo we compute all latent variables #
DO
# E weight (eigen component) Eᵀf. The norm is an eigenvalue #
IF REAL norm = NORM (T e * f);
norm > small real
THEN COLVEC eigen = (T e * f) / norm;
# Compute latent variable t from E factor score #
COLVEC t := (e * eigen) / NORM (e * eigen);
# Use factor loadings p, q to partial out t from e, f (deflation) #
COLVEC p := T e * t, q := T t * f;
e -:= t * T p;
f -:= t * T q;
# Build matrices #
eigens := eigens BEFORE eigen;
new set := new set BEFORE p; # P #
new consts := new consts ABOVE q # Qᵀ #
FI
OD;
Vervolgens berekenen we de vector β. Hier voorspellen we de waarden van de componenten om de nauwkeurigheid van ons PLS-model te controleren. Het inproduct van de bèta-vector met een nieuwe reeks metingen voor één enkel monster levert de voorspelde waarde op. Je zult iets willen doen met de resultaten, hier drukken we ze bijvoorbeeld af.
MATRIX projected set = T new set * eigens; COLVEC beta = eigens * INV projected set * new consts; # Compare to PCR # MATRIX predictions := mean set * beta + col mean consts; print matrix (pls beta, 0); print matrix (predictions, 0)
Als laatste stap presenteer ik een beter doordacht algoritme voor PLS1. De tweede regel in de code hierboven lost in feite een lineair stelsel van vergelijkingen op. In de praktijk zijn er, gelet op de eindige precisie van computers, meer robuuste methoden om een stelsel van lineaire vergelijkingen op te lossen dan inversie van een matrix. Hieronder volgt een gangbare aanpak voor het oplossen van lineaire vergelijkingen, waarbij gebruik wordt gemaakt van singuliere-waarde-decompositie (SVD).
# Solve the linear equation for beta in eigen-space,
which may be more accurate than INV #
OP SOLVE = (MATRIX a, COLVEC b) COLVEC:
BEGIN MATRIX u, v, VECTOR s;
svd decomp (a, u, s, v);
FOR k FROM 2 TO UPB s # Ignore short eigenvectors #
DO IF s[k] / s[1] < 1e-15 # cf. Python numpy package #
THEN s[k] := 0
FI
OD;
CV svd solve (u, s, v, b[:, 1])
END;
PRIO SOLVE = 9;
COLVEC beta := eigens * projected set SOLVE new consts; # Compare to PCR #
MATRIX predictions := mean set * beta + col mean consts
Hiermee is deze demonstratie ten einde. Het is je wellicht opgevallen dat de SVD-, PCR- en PLS-routines daadwerkelijk deel uitmaken van de standaardomgeving van Algol 68 Genie. De broncode van die runtime-routines is gebaseerd op de GNU Scientific Library, die op haar beurt BLAS implementeert. Die code is geschreven in C, dus nu ik Partial Least Squares in niet één maar twee pre-Python-talen heb gepresenteerd, beschouw ik mijn belofte als vervuld.
