Iteracyjne poprawianie rozwiązań LZNK obliczonych algorytmem równań normalnych*

ROCZNIKI POLSKIEGO TOWARZYSTWA MATEMATYCZNEGO Seria III: MATEMATYKA STOSOWANA XXXI (1989)

An d r z e j Ki e ł b a s i ń s k i

Warszawa

Iteracyjne poprawianie rozwiązań LZNK obliczonych algorytmem równań normalnych*

(Praca wpłynęła do Redakcji 1987.11.11)

0. Wprowadzenie. Rozważamy liniowe zadanie najmniejszych kwadratów (LZNK), Ax ~ b, z daną macierzą A e R m,n i wektorem wyrazów wolnych b eRm. Przy założeniu rank (A) = n < m chcemy wyznaczyć wektor x*eR", który minimalizuje normę euklidesową wektora residualnego b - Ax (w pracy tej rozważamy jedynie normę euklidesową wektorów i normę spektralną macierzy).

Jak wiadomo, zachodzą zależności

jx* = A+ b, A+ = (A TA )-ł AT,

^ (J|r*|| < ||b — Ax|| Vx, r* = b -A x * .

Najmniej kosztownym ze znanych algorytmów rozwiązujących to zadanie jest algorytm ^w n ań normalnych Gaussa. Nazwiemy go w skrócie algoryt-

mem RN. Polega on na wyznaczeniu macierzy B i wektora c,

(2) B := A TA, c := A r b,

oraz rozwiązaniu układu równań normalnych »

(3) Bx = c,

np. algorytmem Cholesky’ego, por. [3] cz. 2.

To przejście od danych {A, b) do wyniku {x} jest jednak na ogół niestabilne, tzn. przy realizacji w arytmetyce numerycznej błędy zaokrągleń mogą zniszczyć nadmierną część informacji o rozwiązaniu, zawartej w da- nych.

Okazuje się, że iteracyjne poprawianie pozwala (bez poważnego zwiększe- nia kosztu i bez stosowania arytmetyki wyższej precyzji) odzyskać istotną

* CPBP 01.01.04.

część utraconej informacji, jeśli spełniony jest warunek postaci (4) vKx2 < 1 , K = K(m, n) ~ mn,

gdzie v oznacza względną precyzję arytmetyki zmiennopozycyjnej, wielkość

(5) x = \\A\\\\A + \\

charakteryzuje zaś wrażliwość zadania Ax ~ b.

1. Błędy zaokrągleń w algorytmie RN. Rozważamy tu standardową aryt- metykę zmiennopozycyjną. Stosując analizę pozornych zaburzeń, por. [6] i [3], wykazujemy, że rozwiązanie x, wyznaczone algorytmem RN, spełnia równanie

(6) < (Ar A + dB)x = (A + <5A)Tb, przy czym macierze zaburzeń mają oszacowania

(7) \\SA\\/\\A\\ ś £i = vKl9 \\AB\\/\\A\\2 ^ ^g2=vK 2,

gdzie wielkości K x, K 2 zależą od sposobu realizacji instrukcji (2) oraz rozwiązywania układu (3).

U w aga 1. Przy pełnej (tzn. nierozrzedzonej) macierzy A, zwykłym sumo- waniu iloczynów skalarnych oraz przy zastosowaniu algorytmu Cholesky’ego można przyjąć

K x = mnl/2/2, K 2 = mn/2. ■

U w aga 2. Dla uproszczenia rozważań przyjmijmy, że stosowany jest wariant algorytmu Cholesky’ego, który zastępuje ujemne lub zerowe elemen- ty główne wielkością — v ||B||. Algorytm taki wyznaczy rozwiązanie x, speł- niające (6) i (7) dla każdego regularnego LZNK. ■

Mnożąc (6) obustronnie przez macierz W =(A r A)-1, otrzymujemy (I + WdB)x = x* + W<5ATb,

po odpowiednim zaś uporządkowaniu napiszemy

x —x* = — WdB(x — x*) — WdBx* + W<5Ar b.

Stąd, przy założeniu

(8) rj2 = e2x 2 < 1,

wynika oszacowanie

(9) ||x -x * || ^ (>?2 £* + >h £)/(!->72)>

gdzie

(10) rjx = x 2, £* = ||x*||, P = ||b||/||A||.

Ireracyjnc poprawianie rozwiązań LZNK ... 93

Oszacowanie (9) jest realistyczne. Musimy się więc liczyć z możliwością, że błąd obliczonego rozwiązania x będzie nie mniejszy niż np.

(11) v s/n x 2{£* + /?)/(1 — v x'nx2).

Jak wiadomo, por. [5], błędy reprezentacji danych JA, b] w arytmetyce numerycznej mogą zaburzyć rozwiązanie zadania Ax ^ b co najwyżej o wielkość

(12) v J n x {£* + xg* + /?)/(1 - v v/ n x), gdzie

(13) ^Q* = \\r*\m \\.

Porównując (11) z (12), spostrzegamy, że błąd wytworzony w algorytmie RN będzie prawdopodobnie co najmniej P-krotnie większy, niż przeniesiony błąd reprezentacji danych, gdzie

(14) P = x/( \+xp), n = e*/(P+ł*)-

Wskaźnik x wrażliwości zadania bywa często bardzo duży. Algorytm RN traci więc istotną część informacji, zawartą w danych, gdy iloczyn xp jest mały. Będzie tak np. w przypadku układów słabo sprzecznych (g* <ś /?) oraz układów o wielkich rozwiązaniach (ę* ^ xfi).

Twierdzenie 1. Błędy zaokrągleń w algorytmie RN mogą powodować x razy większą niedokładność wyniku niż błędy reprezentacji danych. Algorytm RN jest więc numerycznie niestabilny w pełnej klasie regularnych LZNK. ■ Iloczyn xp może być jednak rzędu wielkości wskaźnika x, a mianowicie wówczas, gdy równocześnie zachodzą oba następujące przypadki:

— układ jest silnie sprzeczny (g* ~ /?),

— rozwiązanie x* jest małe (<!;* ~ fi).

Gdy występuje koniunkcja obu tych przypadków, wtedy P jest niewielką liczbą, straty dokładności w algorytmie RN są małe. Okazuje się, że w tym szczególnym przypadku algorytm RN jest nie tylko stabilny, ale także numerycznie poprawny. Wynika to z następującego twierdzenia.

Twierdzenie 2. Niech będą spełnione nierówności (15) tj3 = vK 3x2 < 1, K3 = max(Kt , K 2),

(16) n = e*/W + (*)> 0.

Wektor x, spełniający równanie (6) z warunkami (7), spełnia równocześnie równanie

(17) (A + 9A)r [b — (A + 9A) x] = 0,

gdzie SA eRm,n ma oszacowanie

(18) ||8A ||/||A ||«»X 3/[> (l-»l3)].

D ow ód. Równość (6) zapisujemy w postaci

(19) ATr = y,

gdzie

r = b — Ax, y = dBx — <5Ar b.

Z (9) wynika oszacowanie ||x|| ^ + ri3{^* + ^)/(l — rj3), a zatem

||y|| < v x 31|A|| (||A|| ||x|| + ||b||) ^ vK3 ||A|| (||A|| ||^x*|| + ||b||)/(l - r,3) =

= v ^ 31|A|| ||r*||/[/i(l — f/3)].

'Definiując macierz >9A zależnością 0,

- m r 2yT, otrzymujemy oszacowanie, por. (1),

PA || = ||y||/||f|| < ||y||/||r*|| < vX3||A||/D i(l-iy3)], czyli (18). Równocześnie z (19) wynika

(A + 9A)Tr = 0,

gdy r = 0, gdy r # 0,

a więc również

(A + SA)T [b — (A + SA) x] = (A + SA)T [r — 9Ax] = 0. .

Równość (17) oznacza, że x jest rozwiązaniem LZNK (A + $ A )x ~ b . O numerycznej poprawności algorytmu RN możemy jednak mówić tylko wówczas, gdy zagwarantowana jest względna małość macierzy SA.

Musimy w tym celu zastąpić warunki (15), (16) warunkami silniejszymi, np.

tak, jak w następującym twierdzeniu.

Twi er dz en ie 3. Algorytm RN jest numerycznie poprawny z charakterystyką kumulacji K = 20K3, K 3 = m a x ^ , K 2), w podklasie regularnych LZNK określonej warunkami

(20) p = ||r*||/(||b|| + 1|A|| ||x*||) ^ 0.1,

(21) rj3 = vK3x 2 ^ 0.5.

D ow ód wynika natychmiast z twierdzenia 2. Warunek (20) zawiera w sobie koniunkcję warunków:

— układ silnie sprzeczny (g* > >5/10),

— rozwiązanie „małe” (£* < 10^* < 10)3).

Iteracyjne poprawianie rozwiązań LZNK ^... 95

2. Porównanie algorytmu RN z algorytmem Goluba-Householdera. Algo- rytm Goluba-Householdera, w skrócie: algorytm GH, opisany w [5] i [3], jest numerycznie poprawny z charakterystyką kumulacji KCH = 4mn3/2 w pełnej klasie regularnych LZNK. Wyznaczone tym algorytmem rozwiązanie x zadania Ax ~ b spełnia więc równość

(A +

3A)T [(b

SA)

gdzie

||<5A||/||A|| ^ vK gh, ||<5b||/||b|| < vKGH1\ A . Jeśli spełniony jest warunek postaci

(22) <pGH = vKg h x < 1, to zachodzi oszacowanie

llx-x*|| ^ <pGHię*+ xQ *+ py{i-(pGH).

Twierdzenia 1, 2, 3 poprzedniego rozdziału wyjaśniają, że algorytm RN zapewnia podobną jakość rozwiązania x jedynie w podklasie LZNK, okreś- lonej warunkami typu (20), (21). Należy przy tym zauważyć, że warunek (21) jest znacznie silniejszy, niż warunek (22). Zakres sensownej stosowalności algorytmu RN jest więc istotnie węższy, niż zakres stosowalności algorytmu GH.Algorytm RN jest natomiast istotnie tańszy, niż algorytm GH, gdy

stosunek n/m jest mały (opms s koszt 1-go mnożenia i dodawania):

— koszt RN ^ j n 2(m + n/3) opms,

— koszt GH ~ n2(m — n/3) opms.

Drugą zaletą algorytmu RN jest małe obciążenie pamięci komputera.

Jeśli chcemy zachować nietknięte dane początkowe zadania {A, b), to algo- rytm RN wymaga jedynie ~ n2/ 2 miejsc dodatkowych pamięci (by przecho- wać macierz B, a następnie czynnik Cholesky’ego L), natomiast algorytm GH wymaga w tej sytuacji dodatkowych ~ mn miejsc pamięci.

3. Iteracyjne poprawianie rozwiązania LZNK, otrzymanego algorytmem RN. Rozważamy następujący algorytm

(23) B : = Ar A, LLr : = B, xo := 0 , for

i :

(24)(25) (26)

g i : = Ar (b —Ax,), df : = (LLr)_1gi, xi+i := X,- + d,-,

który nazwiemy w skrócie algorytmem RNIP. Druga instrukcja (23) oznacza tu obliczenie czynnika L (macierz trójkątna dolna) rozkładu Cholesky’ego

macierzy B. Instrukcja (25) powinna być wykonana przez rozwiązanie dwu układów trójkątnych

Zauważmy, że obliczony wektor jest identyczny z rozwiązaniem x, obliczonym w algorytmie RN. W następnym rozdziale pokażemy, że następ- ne wyrazy ciągu jxfc), k = 1, 2, mogą być lepszymi (niż xx) przybliżenia- mi rozwiązania x * = A + b. Koszt instrukcji (23) szacuje się wielkością

~-^n2(m + n/3) opms, a koszt obliczenia każdego następnego przybliżenia xi+i -wynosi ^n(2m + n) opms. Jeśli wystarczy wykonanie tylko kilku kro- ków iteracyjnego poprawiania, to koszt całkowity algorytmu RNIP jest tylko nieznacznie większy, niż koszt algorytmu RN. Realizacja RNIP wymaga przeznaczenia jedynie dodatkowych ~ n2/ 2 miejsc pamięci na przechowanie dolnego (lub górnego) trójkąta symetrycznej macierzy B, a następnie czynni- ka trójkątnego L (tak jak w algorytmie RN).

4. Błędy zaokrągleń w algorytmie RNIP. Analiza błędów zaokrągleń w instrukcjach (23)-(26) pozwala zapisać równość

z macierzami zaburzeń <5, A, Ó, A, d, B, które spełniają oszacowania

gdzie wielkości K t, i = 0, 1, 2, zależą od sposobu realizacji instrukcji algoryt- mu.

U w aga 3. Przy założeniach, określonych w uwadze 1, charakterystyki K lt K 2 mają wskazane tam wartości, a K 0 = n312. ■

Poza wielkościami £,, / = 0, 1, 2, określonymi w (28) będziemy posługiwali się w dalszej analizie również wielkościami

(29) <p*=£/X, ^ = e , x 2 (/ = 0, 1, 2).

Dla wygody zapisu oszacowań będziemy wprowadzali dalsze analogiczne wielkości: sh cph rjt dla 1 = 3, 4 ,..., dające się wyrazić przez wielkości pierwotne ze wskaźnikami 1 = 0, 1, 2.

Udowodnimy obecnie

Twie rdz eni e 4 Jeśli jest spełniony warunek Lh = g,, LTd,= h .

(27) (Ar A + Ą B) (xf + i - X,) = (A + A)r [b - (A + A) x,]

(28) ||óiA ||/||A ||^ £o = vK0, ||<5f A||/||A|| ^ = vKlf

\\Ai B \ m \ \ 2 ^ s 2 = vK2,

(30) iy4 = 2rj2 + t]1 +<p0 (l + <Pi) < 1,

Iteracyjne poprawianie rozwiązań LZN K ... 97

to zachodzi nierówność

(31) rj5 = lri2 + rj1 + (p0{\ + ę J ] /{ l-r j2) < t]Ą,

a wyrazy ciągu |xk}, generowanego w algorytmie RNIP, spełniają oszacowania (32) l|xk-x * || <

ykQ*

gdzie

^ lk = rh {l-rjk5)/(l-rj4),

ctk = <p0(l + ^ i) ( l- ^ 5 ) /( l- ^ ) + ^5-

D ow ód. Mnożąc obie strony (27) przez macierz W = (ATA)-1 i porząd- kując odpowiednio składniki, otrzymujemy równość ^

(I + W J{ B) (x,- + i - x*) = WĄ- B (x,- - x*) - A + di A [(xf — x*) + x*] + + W^i Ar [r* — (A + Si A) (Xj — x*) — <5* Ax*].

Wynika stąd oszacowanie

l|Xj

- x*|| <

]5 ||x,- - x*|| + [V ( 1 -

2)] £*.

Uwzględniając zależności ||x0 —x*|| = oraz rjs < 1, otrzymujemy stąd (32) i (33) . .

Z twierdzenia 4 wynika numeryczna stabilność algorytmu RNIP, por.

(12), w klasie LZNK, ograniczonej warunkiem nieco ostrzejszym niż (30).

Przy t/4 ^ 0.5 zachodzi np. oszacowanie

(34) lim sup ||Xj — x*|| ^ v 2max (K0, K x) x [£* + x q*].

i GO

Pokażemy teraz, że przy nieco tylko silniejszym założeniu algorytm RNIP jest numerycznie poprawny. Udowodnimy w tym celu następujące twierdze-

nie.

Twie rdzenie 5. Jeśli spełniony jest warunek

(35) riu = 2rj2 + rjl + {2(pQ + (pl)(l + ril) < 1,

to wyrazy ciągu {xk}, sugerowanego w algorytmie RNIP, spełniają równości (36) (A + Ak A)r [b - (A + Ak A) xk] = 0 (k = 1, 2, ...)

przy czym macierze Ak A są oszacowane nierównościami (37) \\Ak A||/||A|| < ^ e f + e l+ e n ,* ^ "

7(1

gdzie wielkości e8, e9, e10 (określone w (50), (57), (58)) spełniają przy dostate- cznie małym vx2 równości przybliżone

e8 ~ £i» e9 ~ eio ~ el 4‘£2*

7 - Matematyka Stosowana t. 31

(38)

D ow ód. Przyjmijmy k — i+ l i przekształćmy (27) do postaci I (A + <5A)r [b — (A + <5 A) x] = y,

gdzie

x — xk, óA — Sfę _ | A, óA — <5k _ | A,

y = flBfXk-Xfc-!), 0B = B - <5Ar A -(A T + (5A7’)M .

Przedstawimy teraz ideę dalszego dowodu. Równość (38) zastąpimy rów- nością

gdzie macierze AA, AA powstaną odpowiednio z macierzy <5A, <5A przez uzupełnienie ich składnikami, skonstruowanymi z pewnego rozkładu wektora y. Konstrukcja ta stanowi najdłuższą część dowodu. Lemat Stewarta (por.

[4]) wykazuje, iż z (41) wynika istnienie takiej macierzy AA, że spełnione są zależności

Warunek (35) jest silniejszy, niż warunek (30), więc norma ||xk — xfc_ x11 ^

^ ||xk — x*|| + ||xk_ 1 — x*||, ma oszacowanie, wynikające z (32), (33). Tym samym uzyskamy z (40) oszacowania ||y|| oraz konstruowanych zaburzeń, AA i AA. Kładąc Ak A = AA, uzyskamy z (42) zależności (36), (37), co zakończy dowód. ■

Przechodzimy obecnie do konstrukcji i oszacowania zaburzeń AA i AA.

Normę ||xk — xk_ x|| oszacujemy przez wielkości £, ^q, zdefiniowane równoś- ciami

Z (33) i (35) wynika txk < max [ę>0(l + <Pi)/(l — V*), 1] = 1, a zatem napi- szemy

(41) (A-MA)r [b —(A + dA)x] = 0,

(42) (A + AA)T [b — (A + A A) x] = 0,

||JA ||2 «|M A ||2 + |MA|p.

Iteracyjne poprawianie rozwiązań LZNK ... 99

Wykorzystując (32) oraz uzyskane oszacowania q* i £*, otrzymujemy oszaco- wanie

(44) llXfc-Xk-ill < Ae + ££,

przy czym współczynniki A, £ wyrażamy wzorami (por. (33)) A = IrjJil-rjj),

(45) c = 2(p0 (1 -I- 2<Pj)/(1 - rj7) + [1 + rj5 (1 + ćp)/{ 1 - <p)] 7(1 - rfs), gdzie

(46) <£ = ę>o(l + <pi)/(l-»/4), r\n = 2r]2 + rjl + 2ę>0(l+ <Pi)-

(Zauważmy, że warunek (35) gwarantuje, że wielkości rfo, q> i rj4 są mniejsze od 1.) Z (40) i (44) wynika oszacowanie

(47) llylK«6l|A||29,

9

xe + u .

Przedstawmy wektor y w postaci sumy:

y = f+ h ,

(48) , f = 0, h = 0, gdy 5 = 0,

f = yA^/5, h = yc m gdy 5 > 0.

Definiujemy macierz AA równością (49) AAT =T ,^AT —f ||r|r 2rT,

<5 A 7

gdy r ^ 0, gdy r = 0.

Macierz ta spełnia oszacowanie, por. (47),

(50) IMAII/HAII = 6 1[1 + 296/ ( 1 -i?7)], a z (49) wynika zależność

(51) (A + dA)Tr = h.

Macierz A-M A jest regularna, gdyż zachodzą nierówności, por. (35), (52) ||A+ JA|| ^ q>8 = e8 x = q>x [1 + 2?y6/(l - f /7)] < 1.

Zatem wektor h możemy napisać w postaci

(53) h = (A-MA)r h, h = [(A-MA)r] + h.

Uwzględniając (48), (47) i (52), otrzymujemy oszacowanie

(54) l|h|| £61|A|| ^Cx/(1 — (Pg)-

Łącząc (51) i (53) oraz definiując macierz AA równością

<5A, gdy f = 0, (55) AA = <5A + h£ 2 x , gdy £ > 0, otrzymujemy równość (41), a z (54) i (45) oszacowanie

(56) |MA||/||A|| +

gdzie

(57) «9 = e„[l + 2^7(l + 2<p i)/[(1 — J/4.) (1 — <Ps)]] >

(58) £,o =E6[l + IJ5(l + cp)/(l-<p)]/(l-<p8).

Oszacowania (50), (56) oraz nierówność (42) pozwalają napisać oszacowanie (37). Z (40), (50), (57), (58) wynikają podane w twierdzeniu 5 charakterystyki asymptotyczne dla e8, eg, el0. ■

Skomplikowana zależność wielkości r\u przy l > 2 od wielkości pier- wotnych zaciemnia sens otrzymanych charakterystyk. Zapisujemy więc nastę- pującą, osłabioną, ale bardziej przejrzystą wersję twierdzenia 5.

Tw ie r dz e nie 6. Jeśli dla 1 = 0, 1, 2 spełnione są warunki (59) vKtx = <Pi < ę < 0.03, vKtx 2 = y]x ^ rj ^ 0.2, to zachodzą nierówności

(60) fj = rj5 ^ (2rj+ I.03ę>)/(1 — rj) < 2.5rj + l.3ę < 0.54,

a wyrazy ciągu {xk}, generowanego w algorytmie RNIP, spełniają równania (A + Ak A)T [b — (A + Ak A) xfc] = 0 (k = 1, 2, ...),

przy czym macierze Ak A mają oszacowania

|M»A||/||A|K v {s /Kl + Kt [l+(8^ + 6V) /( l- 3 f))] +

+ (K0 + A1+A:2)tl+ (3 l/ + 2<)>)/(l-);)]^-‘- 1/(l-^ )} . . Twierdzenie 6 oznacza, że algorytm RNIP jest numerycznie poprawny z charakterystyką kumulacji K = 5.5 yjK^ + K f w klasie regularnych LZNK, ograniczonej warunkiem (59).

5. Wnioski końcowe. Proces RNIP może być efektywnie kontrolowany.

Z dowodu twierdzenia 5 wynika, że maksymalna graniczna dokładność jest osiągana w RNIP, gdy norma poprawki, ||dfc|| % H x ^-x ^ ll, jest rzędu

wielkości

2vx(Kl x ||b —Axk||/||A|| + K 0 ||xk||).

Iteracyjne poprawianie rozwiązań LZNK ... 101

Jeśli więc potrafimy wyznaczyć (choćby w przybliżeniu) wielkości

= 2vxmax(K0, K t), 02 = x/||A|| = ||A + ||,

to wystarczy w algorytmie RNIP zastąpić instrukcję (24) instrukcjami (61) je, := b —Ax,-, w:=||r,.||, z := ||x f||, := ||d£_ 1||

I if d ^ 0i *(O2*w-\-z) then go to END else g ,:= A Tr,-.

Sposób prostszy i zapewne równie skuteczny polega na uzupełnieniu instrukcji (23) instrukcją d d : = 0, a instrukcji (26) instrukcjami

(62) d := ||d,||, if^ i :> 0 Ad > dd then go to END else dd := d/2.

(Przestajemy poprawiać, gdy nowa poprawka nie jest wyraźnie mniejsza od poprzedniej.)

Eksperymenty numeryczne potwierdzają opis teoretyczny procesu, p. [2].

Proces RNIP ma oczywiście inny, węższy zakres zastosowań, niż algo- rytm Bjórcka [1]. Może służyć tylko do odzyskania części informacji, utraconej w algorytmie RN, gdy spełniony jest warunek typu (4).

Prace cytowane

[1] A. Bjfirek, Iterative refinement of linear least squares solutions I, BIT 7 (1976), 257-278.

[2] A. G ołąb , Proces RNIP, Praca magist., Wydz. MiM UW, 1987.

[3] J. i M. Jan k ow scy, Przegląd metod i algorytmów numerycznych, cz. .1, WNT, Warszawa 1981.

M. D ryja, J. i M. Jan k ow scy, Przegląd metod i algorytmów numerycznych, cz. 2, WNT, Warszawa 1982.

[4] G. W. Stew art, On the perturbation on pseudo-inverses, projections and linear least squares problems, SIAM Rev. 19 (1977), 634-662.

[5] J. Stoer, Wstęp do metod numerycznych, t. 1, PWN, Warszawa 1979.

[6] J. H. W ilk in son , Błędy zaokrągleń w procesach algebraicznych, PWN, Warszawa 1967.