Porównanie algorytmów minimalizacji w zastosowaniu do analizy złożonych sygnałów chromatograficznych

W itold ILEW IC Z Politechnika Śląska

P O R Ó W N A N IE A L G O R Y T M Ó W M IN IM A L IZ A C JI W Z A S T O S O W A N IU D O A N A L IZ Y Z Ł O Ż O N Y C H S Y G N A Ł Ó W C H R O M A T O G R A F IC Z N Y C H

S treszczenie. W artykule porów nano 3 iteracyjne algorytm y m inim alizacji, zaim plem entow ane w środow isku M atlab, w zastosow aniu do analizy złożonych sygnałów chrom atograficznych m etod ą najm niejszych kw adratów : Gaussa- N ew tona, Levenberga-M arquardta oraz m etodą sim plex N edlera-M eada. Podsta­

w ow e kryteria porów nania były następujące: czułość algorytm u na dokładność wstępnej oceny param etrów m odelu, szybkość i zbieżność algorytm u oraz przy­

datność algorytm u do autom atycznej analizy sygnałów chrom atograficznych.

Badania przeprow adzono m etodam i sym ulacyjnym i dla typow ych m odeli stoso­

w anych w chrom atografii.

C O M P A R IS O N O F M IN IM IZ A T IO N A L G O R IT H M S F O R P U R P O S E O F A N A L Y S IS O F C O M P L E X C H R O M A T O G R A P H IC SIG N A LS

S u m m a ry . In article there are com pared 3 iterative algorithm s o f m inim ization included in M atlab environm ent, for purpose o f analysis o f com plex chrom atographic signals. These are: G auss-N ew ton m ethod, Levenberg- M arquardt m ethod and N edler-M ead Sim plex m ethod. Basic critcrions o f com parison w ere sensitivity o f m ethod to starting param eter vector, speed, convergence o f algorithm s to proper solution and usefulness o f algorithm to autom atic analysis o f chromatogram s.

1. W stęp

W chrom atografii uzyskiw ane sygnały m ają postać sum y pików i nazyw a się je chrom atogram am i. C hrom atogram je s t to cyfrow y zapis przebiegu sygnału z detektora chrom atografu w funkcji czasu, dyskretny w czasie (próbkowanie) oraz w w artościach (kwantowanie). Podstaw ow ym i elem entam i sygnału chrom atograficznego są piki, z których każdy odpow iada pew nem u składnikowi w m ieszaninie rozdzielanej w chrom atografie. Gdy rozdzielczość m etody analitycznej je st mała, uzyskiw any sygnał składa się z pików nakładających się na siebie i taki sygnał je s t traktow any w pracy jak o złożony. Celem analizy takiego sygnału jest w yznaczenie położenia i pow ierzchni każdego z pików składowych w sygnale. Pow ierzchnia piku stanowi

podstaw ę analizy ilościowej i pozw ala określić zaw artość danego składnika w analizowanej m ieszaninie, położenie piku pozw ala zidentyfikow ać dany składnik.

P odstaw ow ą m etodą analizy złożonych sygnałów chrom atograficznych je st regresja nieliniow a [1][2], Dla m inim alizow anego w skaźnika postaci:

J(b) = £ ( y ; - )) , (1)

/=!

gdzie: F(b,l,-) stanowi m odel sygnału w postaci sum y pików zależny od w ektora param etrów b, yr}'(U ) reprezentuje w artości próbek zarejestrow anego chrom atograinu, poszukuje się takiego w ektora param etrów b , że zachodzi:

J (b) = m in (J(b )). (2)

N a podstaw ie w ektora param etrów b = [b ,, b ,, .. b J T (q stanow i liczbę pików składowych w sygnale) m inim alizującego w skaźnik ./(b) oraz m odelu pojedynczego piku odtwarza się kształty pików składow ych i wyznacza się ich pow ierzchnię i położenie.

Zadanie estymacji param etrów m odelu realizuje się m etodam i iteracyjnym i ze względu na nieliniow ość m odelu w zględem param etrów . Podstaw ow e m etody stosowane w takich przypadkach to: m etoda G aussa-N ew tona (GN), m etoda Levenbcrga-M arquardta (LM ) oraz m etoda Sim plex N edlera-M eada (SPX), szeroko om aw iane w literaturze (np. [3],[4]). Dwie pierw sze są m etodam i gradientow ym i, przy czym m etoda G aussa-N ew tona opiera się na kwadratowej aproksym acji w skaźnika jakości, natom iast m etoda Levenberga-M arquardta je st połączeniem m etody naj­

w iększego spadku oraz m etody G aussa-Newtona. Trzeci z rozw ażanych algorytm ów - Sim plex - je st algorytm em bezgradientow ym , w którym odpow iednio m odyfikuje się w ielościan o «+1 w ierzchołkach w przestrzeni param etrów , zw any sym pleksem .

W literaturze opisuje się w ady i zalety w ym ienionych algorytm ów . A lgorytm G aussa-N ew tona je st szybko zbieżny, jednakże je st czuły na dokładność doboru startow ego w ektora param etrów m odelu. Jego rozw inięciem je st algorytm L-M , w którym w początkow ych iteracjach stosuje się algorytm najw iększego spadku, który w pobliżu m inim um zastępuje się algorytm em G-N. A lgorytm ten je st jednym z najszerzej stosowanych w różnych problem ach optym alizacji. A lgorytm „Sym plex”

m a trzy podkreślane w literaturze zalety w stosunku do m etod dyskutow anych wcześniej. Po pierw sze - je s t prosty - nie w ym aga w yznaczania pochodnych w skaźnika jakości w zględem param etrów m odelu (gradientu i hesjanu), po drugie - zbieżność rzadko następuje w m inim um lokalnym , po trzecie - algorytm m oże być stosowany do m inim alizacji nieciągłych funkcji.

Ogólne własności w ym ienionych algorytm ów m inim alizacji są znane, jednakże w arto spojrzeć na problem stosow ania tych algorytm ów w kontekście autom atyzacji procesu analizy złożonych sygnałów chrom atograficznych. K onieczność elim inacji człow ieka z procesu analizy sygnału w ynika z potrzeby analizow ania coraz bardziej skom plikow anych m ieszanin, często prow adzących do sygnałów o liczbie pików rzędu kilkuset, z których w iele nakłada się na siebie, zw łaszcza w dziedzinie analiz w ochronie środowiska [5] (np. badania czystości pow ietrza). Interaktyw na analiza takich sygnałów z udziałem człow ieka je s t niecelow a, w obec czego pow staje problem

w iarygodności w yników otrzym yw anych z zastosow aniem algorytm ów autom atycznej analizy, której jed n y m z etapów je st ocena param etrów nakładających się pików m etodą regresji nieliniow ej. Stąd wyniki pracy pozw alają w nioskow ać o przydatności ww. algorytm ów do autom atycznej analizy złożonych sygnałów chrom atograficznych.

2. T esto w an ie alg o ry tm ó w G N , L M i SPX

Testy algorytm ów przeprow adzono w środow isku M atlab dla 9 sym ulow anych chrom atogram ów z dwom a nałożonym i pikam i składowym i w g m odelu gausso­

wskiego:

2w2 ^{(3 )}

oraz 9 sygnałów z dw om a pikam i asym etrycznym i w g m odelu Frasera-Suzkiego:

f ( h , t⁰,w ,a ,t) = h-cxp ■ ln(2)

ln 1 +

w- a ^{(4 )}

gdzie: h - w ysokość piku, w - szerokość piku, t0 - położenie piku, a - w spółczynnik asym etrii piku. W porów naniu do piku gaussow skiego, m odel w g Frasera-Suzukiego charakteryzuje się dodatkow ym param etrem a. Zm iany param etru a w pływ ają na zm ianę stopnia asym etrii piku i tak dla a > 0 pik staje się asym etryczny w prawo, dla a < 0 pik staje się asym etryczny w lewo, przy czym stopień asym etrii jest proporcjonalny do a . W ygenerow ane sygnały przedstaw iono na rysunkach la ) i lb).

Procedura testow ania algorytm ów obejm ow ała 100-krotną estym ację param etrów każdego z sygnałów testow ych za pom ocą każdego z 3 algorytm ów, przy czym za każdym razem sygnał testow y zaszum iano szum em gaussow skim N(0,0.05), a startowy w ektor param etrów za każdym razem był generow any na podstaw ie praw dziw ego poprzez zakłócenie go na poziom ie m aksym alnym 50% wartości param etrów zm iennym i losowym i o rozkładzie prostokątnym . Szacuje się, że w praktyce dokładność wstępnej oceny w ektora param etrów m odelu dla sygnałów chrom atograficznych jest na poziom ie 10-20%, jed n ak przyjęty tu 50% próg zakłóceń startowego w ektora param etrów m iał uw ypuklić ew entualne różnice w w ynikach uzyskiw anych dla poszczególnych algorytm ów .

Badania algorytm ów przeprow adzono w środow isku M atlab, w ykorzystując go­

tow e funkcje realizujące ww. algorytm y. S ą to: funkcja ‘n linfit’ ze Statistics Toolbox, w której zaim plem entow ano algorytm G aussa-N ew tona, funkcja ‘lsqcurvefit’

z O ptim ization Toolbox z zaim plem entow anym algorytm em Levenberga-M arquardta, funkcja ‘fm insearch’, za pom ocą której realizuje się algorytm Sim plex N edlera- M eada. Funkcje te realizują m inim alizację sum y kw adratów odchyleń m odelu od sygnału bez ograniczeń, to znaczy w kolejnych iteracjach nie nakłada się ograniczeń na m ożliw e w artości param etrów m odelu.

3. W y n ik i

W yniki dla sygnałów z gaussow skim m odelem pików przedstaw iono na rysunkach 2 i 3. N a rysunku 2 zaprezentow ano wyniki 900 realizacji każdego z testow anych algorytm ów na sygnałach testow ych z rysunku la.

a)

[1] A tAv=1.3 h ^ . ^ 1

>

A A

/ A ‘

y_ J \ __

[4] 4 t * = 1 . 3 h , * j = 2

f \A / U

V \

[7J A tAv=1.3 ł ^ /h ^ lO

A

.< i

i \

[2] A tAv=0.93 h1/h2=1

M / \ / / \ ■

[5] A t/w =0.93 h,/h2=2

/«.A

[3] A tAv=0.47 h,/h 2=1

\

/ A

' A \

18) A t*i= 0 .9 3 h,/h2=10

t \

A

|9] A łAv=0.47 h1/h2=10

A

/ \

b)

[1] 8 ^ 0 . 2 8 2 = 0 . 3 . ^ ^ 2 = 1 [2 ] a 1= 0 .4 a 2 = 0 .5 . h t /h 2= 1 [3 ] a ^ O . 6 8 ^ 0 . 7 . h .,/h 2= 1

A A

/ Y \

J JVA__

[4] 8 ^ 0 .2 8 2 = 0 . 3 . h1/h 2= 2

A

;V\

[7 ] a 1s 0 . 2 a 2= 0 .3 . h 1/h 2 * l 0

i 1/ \ i ]

!

i \

/ \

f \

i \ J X

i \

R y s.l. a) Sygnały testowe od [1] do [9] w postaci dw óch nałożonych pików wg m odelu gaussow skiego; At/w - stosunek różnicy położeń pików składow ych do szerokości piku, h]/h2 - stosunek w ysokości pików składow ych; b) sygnały testow e od [1] do [9] w postaci dw óch nałożonych pików wg m odelu Frasera- Suzukiego; hi/h2 - stosunek w ysokości pików składow ych, ai, a2 - w spółczynniki asym etrii dla piku pierw szego i drugiego; kolor niebieski - sygnał; kolor czerw ony - piki składow e

Odcięta każdego punktu w idocznego na w ykresach na rysunku 2, oznacza odległość punktu startowego w przestrzeni param etrów od punktu optym alnego, w yrażona jak o procent norm y punktu optym alnego (w przestrzeni param etrów).

R zędna każdego punktu w yraża odległość uzyskanego z zastosow aniem danego algorytm u punktu reprezentującego rozw iązanie od punktu optym alnego, także w yrażoną jak o % norm y punktu optymalnego.

GN LM SPX

Rys. 2. Porów nanie zbieżności m etod GN, LM i SPX dla m odelu gaussowskiego Ze w zględu na d u żą zm ienność na osi rzędnych przyjęto skalę logarytmiczną.

Linie w idoczne na rysunku 2 rozdzielają zbiory punktów uzyskane w wyniku sym ulacji na 2 podzbiory. Punkty leżące powyżej tej linii oznaczają przypadki, w których, po w ykonaniu algorytm u, odległość punktu otrzym anego jak o rozwiązanie od punktu optym alnego je st w iększa od odległości punktu startowego do optym alnego przed m inim alizacją. A więc są to przypadki, gdy algorytm okazał się rozbieżny lub uzyskano zbieżność do lokalnego m inim um . Liczba przypadków braku poprawy dokładności rozw iązania w porów naniu do startow ego w ektora param etrów m ożna oszacow ać na ok. 40% w przypadku algorytm u GN, 10% w przypadku algorytm u LM, i 20% w przypadku algorytm u SPX.

Zbiorcze rozkłady czasów w ykonania algorytm ów podczas analizy sygnałów z pikam i gaussow skim i z rysunku la) przedstaw iono na rysunku 3. N ajszybszą m etodą je st m etoda GN. Średni czas m inim alizacji w skaźnika jakości do spełnienia kryterium stopu w ynosi ok. 0.1 ms. W idać jednakże, że w ystępują odstające przypadki w ydłużenia tego czasu, naw et pięciokrotnie. Przypadki te oznaczają najpraw do­

podobniej w ystąpienie niezbieżności algorytm u i osiągnięcie m aksym alnej zadanej liczby iteracji, co kończy realizację algorytmu.

W porów naniu do algorytm ów LM i SPX rozrzut czasu analiz dla algorytm u GN je st kilkakrotnie m niejszy. A lgorytm LM charakteryzuje się lepszą średnią szybkością n iż algorytm SPX, jednakże występuje tu najw iększe rozproszenie w yników . Spow olnienie w stosunku do m etody GN spowodowane jest przełączaniem

algorytm u na m etodę najw iększego spadku, zapew niającą lepszą zbieżność algorytm u LM w stosunku do GN. K osztem spow olnienia uzyskuje się w ięk szą odporność na dokładność doboru wartości startowego w ektora param etrów . Średnio najgorszy je st algorytm SPX, przy czym różnica w szybkości dla m etody SPX i LM je st dużo mniej istotna niż różnice m iędzy algorytm am i SPX i LM a algorytm em GN.

I

0.5

Rys. 3. Porów nanie czasów w ykonania algorytm ów GN, LM i SPX dla m odelu gaussowskiego

W yniki analogiczne ja k dla sygnałów z pikam i gaussow skim i uzyskane dla sygnałów z pikam i wg m odelu Frasera-Suzkiego przedstaw iono na rysunkach 4 i 5.

6 SPX 7

6

5

4

3 ■ 2 • t - 0 -1

-2

-3

R ys. 4. Porów nanie zbieżności m etod GN, LM i SPX dla m odelu Frasera-Suzukiego

0 8

0.7

0.6

0 5

E

“ 0 4 -

0.3 •

0.2

0

ON LM SPX

Metoda

Rys. 5. P orów nanie czasów w ykonania algorytm ów GN, LM i SPX dla m odelu Frasera-Suzukiego

4. P o d su m o w an ie i w nioski

W św ietle przeprow adzonych badań sym ulacyjnych najlepsze wyniki uzyskano dla algorytm u Levenberga-M arquardta (LM). A lgorytm ten okazał się najodporniejszy na niedokładności startow ego w ektora param etrów , dawał najm niejszą liczbę przypadków rozbieżności w porów naniu z pozostałym i dwom a. C zas rozw iązyw ania zadania regresji je st dla algorytm u LM średnio większy w porów naniu z algorytm em GN, lecz m niejszy w porów naniu z algorytm em SPX.

A lgorytm GN zapew niał najkrótszy czas osiągnięcia rozw iązania, natom iast okazał się najm niej odporny na dokładność określenia startow ego w ektora param etrów - średnio najczęściej, w porów naniu z pozostałym i algorytm am i, prow adziło to do przypadków rozbieżności algorytm u lub zbieżności do m inim um lokalnego znacznie oddalonego od punktu optym alnego.

A lgorytm SPX okazał się najw olniejszy spośród 3 badanych algorytm ów. Dla w iększych w ym iarów zadania czas uzyskania rozw iązania m oże być w ielokrotnie w iększy niż dla pozostałych dwóch algorytm ów. W literaturze podkreśla się za to jeg o odporność na zakłócenia, co jednak, w przypadku złożonych sygnałów chrom ato­

graficznych nie potw ierdziło się.

W ydaje się wobec pow yższego, że algorytm LM je s t najlepiej dostosow any do przeprow adzania autom atycznej analizy złożonych sygnałów chrom atograficznych spośród badanych w pracy, co jednocześnie nie oznacza, że je st idealny. Stąd warta rozw ażenia je st kom binacja algorytm u LM z innymi algorytm am i estym acji modeli param etrycznych, np. genetycznym i czy też EM (Expectation-M axim ization).

L1TERATURA

1. Reh E.: Peak-shape analysis for unresolved peaks in chrom atography:

com parision o f algorithm s. Trends In Analytical C hem istry, v o l.14, n o .l, 1995.

2. Scheeren P.J.H., B am a P., Smit H.C.: A Software Package for the Evaluation o f Peak Param eters in an A nalytical Signal Based on a N on-L inear R egression M ethod, A nalytica C him ica Acta, 167,1985, p. 65-80.

3. M arquardt D.W .: A n algorithm for least-squares estim ation o f nonlinear param eters. J. Soc. Ind. Appl. M ath. 11, 1963, p. 431-444.

4. Caceci M .S., Cacheris W .P.: Fitting curves to the data. The sim plex algorithm is the answer. Byte 9, 1984, p. 340-362;

5. Steffen B. i inni: A new m athem atical procedure to evaluate peaks in com plex chrom atogram s. Journal o f C hrom atigraphy A, 1071, 2005, p. 239-246.

Recenzent: D r hab. inz. A dam K ow alew ski, prof. A G H

A b stra c t

In article there are com pared 3 iterative algorithm s o f m inim ization included in M atlab environm ent, for purpose o f analysis o f com plex chrom atographic signals.

These are: G auss-N ew ton m ethod, Levenberg-M arquardt m ethod and N edler-M ead Sim plex m ethod. Basic critcrions o f com parison w ere sensitivity o f algorithm s to starting param eter vector, speed, convergence to proper solution and usefulness to autom atic analysis o f chrom atogram s. C om parison was done during sim ulations, where there w ere analyzed signals w ith two strongly overlapping peaks for gaussian peak model and Fraser-Suzuki peak model. A ccording to results o f sim ulations an algorithm o f Levenberg-M arquardt seem s to be regarded the m ost effective and best suited to purpose o f autom atic analysis o f com plex chrom atographic signals. This algorithm was the m ost robust for inaccuracy in starting param eters vector com paring to two others algorithm s. It w a sn ’t the fastest one, but w asn ’t the slow est one too and in m any cases its speed was com parable to the fastest algorithm - Gauss-Newton. It seems to be prom ising com bination o f LM algorithm w ith other algorithm for param eter estim ation problem s like genetic algorithm s or EM (Expectation- M axim ization).