 |
 |
 |
|
{Plus ça change …}
|
Toen ik in de jaren tachtig studeerde, stond de chemometrie nog in de kinderschoenen. Voor ons, studenten, leek het een esoterisch vakgebied dat alleen in nicheprojecten op afdelingen analytische chemie aan bod kwam. Over het algemeen was de rekencapaciteit van computers nog niet toereikend om de grote datasets te verwerken die we tegenwoordig zien. Chemometrie kon zijn volledige potentieel nog niet laten zien – achteraf gezien was het gewoon weer een vakgebied waar de theorie ver voorliep op de toepassing. Voor serieuze number crunching koos je in die tijd voor een theoretisch vakgebied als quantumchemie of statistische fysica. Dat is tenminste wat ik deed.
Men kon er vroeger alleen maar van dromen dat de rekenkracht zo zou toenemen als in de afgelopen veertig jaar het geval is geweest. We kunnen eindelijk grote datasets verwerken die voorheen onmogelijk te hanteren leken, en hebben nieuwe termen bedacht zoals data science of big data. Chemometrie leunt sterk op multivariate statistiek, wat tegenwoordig wordt beschouwd als een machine learning-techniek. Begonnen in de statistische chemische fysica ben ik uiteindelijk in de chemometrie terechtgekomen en de kans is groot dat ik daar met plezier blijf. Chemometrie is een uitstekende keuze voor mensen zoals ik, die graag werken aan projecten waar chemie, fysica, wiskunde en informatica samenkomen. De afgelopen jaren is er een nieuwe overkoepelende term voor die projecten ontstaan: cheminformatica. Die term bevalt me.
Chemometrische algoritmen worden ondersteund door de vele tools voor de ontwikkeling van machine learning-toepassingen. Ik programmeer de meeste van mijn algoritmen zelf, uiteraard in Python, dat een soort moderne lingua franca lijkt te zijn. En sinds mijn dochter, die nu studeert, mij Jupyter Notebooks heeft uitgelegd, voel ik me au courant en zelfverzekerd. Merk op dat "modern" een relatief begrip is – is Python zelf niet al in de dertig, een millennial om het zo te zeggen? Toch vraag ik me toch af hoe de dagelijkse praktijk er voor die pioniers in de jaren tachtig uit moet hebben gezien.
Algoritmen voor machine learning zijn vaak gebaseerd op toegepaste lineaire algebra. Al veertig jaar geleden waren er volwaardige bibliotheken voor numerieke lineaire algebra beschikbaar, zoals BLAS, LINPACK of EISPACK. De laatste twee zijn in de loop der tijd opgenomen in LAPACK, een referentie-implementatie voor numerieke lineaire algebra. Moderne Python-bibliotheken maken onder de motorkap gebruik van BLAS en LAPACK – denk bijvoorbeeld aan numpy.linalg of scipy.linalg.lapack.
Aangezien LAPACK en zijn voorgangers zijn geprogrammeerd in Fortran, vermoed ik dat veel van het chemometrische werk destijds is geprogrammeerd in oudere versies van die taal. Hier wil ik laten zien hoe sommige van die algoritmen er in ouderwets Fortran uit zouden kunnen hebben gezien. Aangezien het slechts een illustratie is, zijn de algoritmen vereenvoudigd. De code in dit bericht kan worden uitgevoerd met VIF, mijn eigen experimentele vintage Fortran-compiler, waarschijnlijk de enige implementatie die REAL*32 en COMPLEX*64 standaard ondersteunt.
De Moore-Penrose-inverse of pseudo-inverse, wordt vaak gebruikt in de chemometrie, omdat deze een kleinste-kwadraten oplossing kan opleveren voor een onderbepaald stelsel lineaire vergelijkingen. Een dataset die bijvoorbeeld resultaten bevat van minder laboratoriummonsters dan er waarnemingen per monster zijn, is onderbepaald. Gegeven een matrix A, kan de pseudo-inverse A+ eenvoudig worden berekend aan de hand van de singuliere-waarde-decompositie van A.
De singuliere-waarde-decompositie stelt dat elke m×n-matrix A een product is A = U Σ VT, waarbij U UT = 1 en V VT = 1. Dan geldt A+ = V Σ+ UT. Σ+ is een diagonale matrix waarvan de elementen de reciproke, van nul verschillende diagonale elementen zijn van Σ, een diagonale matrix van singuliere waarden. In het algemeen geldt dat A A+ A = A. Wanneer A vierkant en inverteerbaar is, geldt A+ = A-1 en A A+ = 1.
De onderstaande code berekent de pseudo-inverse van een matrix. De subroutine roept de LAPACK-functie DGESDD aan, een singuliere-waarde decompositie volgens een "divide-and-conquer"-methode. Merk op dat de Python-functie numpy.linalg.pinv de functie DGESDD aanroept. BLAS beschikt niet over een routine voor de vermenigvuldiging van diagonale matrices, dus die is in deze routine zelf opgenomen.
subroutine pinv(A, m, n, A inv, U, V T, S, tol, info)
! Compute pseudo-inverse of 'A'.
! 'A' is a 'm'x'n' matrix which is overwritten.
! 'A inv' will contain the 'n'x'm' pseudo-inverse.
! 'U' will contain U, and 'V T' will contain V^T.
! 'S' will contain reciprocal singular values.
! 'tol' is the threshold ratio relative to the largest singular value.
! 'info' is a status indicator. Zero on exit is good.
implicit none
integer i, j, k, m, n, lwork, info, iwork, lwmax
real*8 A(m, n), A inv(n, m), U(m, m), V T(n, n), S(n), tol
! Large arrays are shared between subroutines.
parameter (lwmax = 100000)
common /share/ work(lwmax), iwork(lwmax)
real*8 work, smax
integer*4 iwork
! We sum in quad precision.
real*16 sum
! Work size should be at least 8 * min(m, n).
if (lwmax < 8 * min(m, n)) then
info = 1
return
end if
! Find optimal workspace.
lwork = -1
call dgesdd('a', m, n, A, m, S, U, m, V T, n, work, lwork, iwork, info)
if (lwmax < work(1)) then
info = 2
return
end if
lwork = min(lwmax, int (work(1)))
! Compute SVD.
call dgesdd('a', m, n, A, m, S, U, m, V T, n, work, lwork, iwork, info)
if (info > 0) then
info = 3
return
end if
! Construct S^+.
! Too small singular values cause instability, so reject.
smax = S(1)
do i = 1, n
if (S(i) > tol * smax) then
S(i) = 1 / S(i)
else
S(i) = 0
end if
end do
! Compute pseudo inverse V * S^+ * U^T
! BLAS has no routines for multiplying diagonal matrices.
do i = 1, n
do j = 1, m
sum = 0
do k = 1, m
sum = sum + V T(k, i) * S(k) * U(j, k)
end do
A inv(i, j) = dble (sum)
end do
end do
info = 0
return
end
Ik zou zeggen dat de bovenstaande code zelfs voor iemand die nog nooit Fortran heeft gezien, niet vreemd overkomt. Let wel: in het oude Fortran kan geheugen niet dynamisch worden gealloceerd. Daarom moeten alle arrays in de aanroepende functie worden gedeclareerd en aan de routine worden doorgegeven.
Laten we een eenvoudige toepassing maken met de pseudo-inverse. Ordinary Least Square (OLS-) regressie is het oplossen van A x = y voor x, wanneer er geen exacte oplossing bestaat. In de chemometrische praktijk impliceert experimentele meetfout doorgaans het ontbreken van een exacte oplossing. Als een dataset in de vorm van een m×n-matrix A m proefmonsters zou vertegenwoordigen met n waarnemingen per monster, en y de vereiste uitkomst voor elk monster zou zijn, dan zou x een kalibratiefactor zijn. Dit is een eenvoudig voorbeeld van het "trainen" van een machine-learning "model". Het inproduct van een onbekende reeks van n waarnemingen met x zou de uitkomst voor dat monster in het betreffende experiment voorspellen.
OLS-regressie kan eenvoudig in Fortran worden geprogrammeerd, zoals hieronder wordt getoond. Het berekenen van de pseudo-inverse is voor OLS niet nodig, maar houd er rekening mee dat dit slechts een illustratie is. De subroutine roept de BLAS-functie DGEMV aan, een gegeneraliseerd matrix-vectorproduct.
subroutine ols(A, m, n, x, y, A inv, U, V T, S, tol, info)
! Solve 'A' 'x' = 'y' by computing the pseudo-inverse of 'A'.
! 'A' is a 'm'x'n' matrix which is overwritten.
! 'Ainv' will contain the 'n'x'm' pseudo-inverse.
! 'U' will contain U, and 'VT' will contain V^T.
! 'S' will contain reciprocal singular values.
! 'tol' is the threshold ratio relative to the largest singular value.
! 'info' is a status indicator. Zero on exit is good.
implicit none
integer m, n, info
real*8 A(m, n), A inv(n, m), x(m), y(m), U(m, m), V T(n, n), S(n), tol
! Solve A x = b.
call pinv (A, m, n, A inv, U, V T, S, tol, info)
if (info > 0) return
call dgemv ('n', n, m, 1.0d0, A inv, n, y, 1, 0.0d0, x, 1)
return
end
Aangezien chemici vaak datasets produceren die veel correlaties bevatten, is het een betere aanpak dan OLS om de regressie te baseren op dominante eigenvectoren van de covariantiematrix van je dataset. Lange eigenvectoren duiden op richtingen met grote variatie in een dataset, korte eigenvectoren duiden op richtingen met grote correlatie. De beste voorspelling kan worden vastgesteld door de voorspellingsfout te minimaliseren als functie van het toenemende aantal in aanmerking genomen eigenvectoren. De voorspellingsfout neemt doorgaans eerst af (vanwege aanvankelijke underfitting), maar stijgt na een minimum weer wanneer correlaties de kwaliteit van de voorspelling steeds meer gaan beïnvloeden. Dit proces wordt om voor de hand liggende redenen dimensiereductie genoemd, en het op deze manier opstellen van een regressie wordt Principal Component Regression (PCR) genoemd.
Je zult opmerken dat OLS in feite PCR is waarbij alle eigenvectoren zijn meegenomen. Daarom is OLS in de chemometrie geen a priori goede benadering. Als je het nog beter wilt doen dan met PCR, kun je PLS – Partial Least Squares-regressie – overwegen. De vervolgvraag is dan hoe PLS eruit zou zien in een taal van vóór Python? Dit zal het onderwerp zijn van een toekomstige post.
