Analiza numeryczna typowych zadań z unitarnym przekształceniem

ANDRZEJ KIELBASIŃSKI (Warszawa) i ^KRYSTYNA ZIĘTAK (Wrocław)

Analiza numeryczna typowych ^zadań z unitarnym przekształceniem

Householdera

(Praca ^przyjęta do druku 24.4.1975)

O. ^Wstęp. W pracy tej przedstawiamy dowód numerycznej poprawności unitar- nego przekształcenia Householdera. Rozważamy różne arytmetyki numeryczne rzeczywiste lub zespolone oraz typowe !Ub ogólne warunki definiujące transformację.

Pewne szczegóły tej analizy wydają się nam nowe, zasadnicza jednak jej struktura pochodzi z pracy J. H. Wilkinsona [6].

Główną zatem ^przyczyną naszej publikacji jest chęć spopularyzowania tej ważnej

transformacji oraz metod analizy Wilkinsonowskiej.

Numeracja wzorów jest oddzielna w każdym z paragrafów. Cytując wzory poda- ne w innym paragrafie podajemy dwie liczby: numer paragrafu i numer wzoru.

1. Zespolone i rzeczywiste przekształcenia Househołdera. Niech Kn oznacza zespo-

loną lub rzeczywistą kartezjańską przestrzeń wektorów, utożsamianych z macierzami jednokolumnowymi, z normą euklidesową:

Unitarnym przekształceniem Househo/dera nazywamy przekształcenie dające się zapisać w postaci

(I) P = 1-w.2-w 8 , llwll ^{= 1.}

Celowo zastosowaliśmy tu zapis postaci w . ^{2 .} wH' ^gdyż jest to poprawnie ^określony iloczyn macierzy wymiarów n x I, I x 1, I x n (przez BH rozumiemy macierz ^sprzę­

żoną z ^macierzą B).

Przekształcenie (I) przyporządkowuje wektorowi x ^{E Kn} jego odbicie zwierciadla- ne y względem hiperpłaszczyzny prostopadłej do wektora w. ^Istotnie,

-y ⁼ Px ⁼ x-r.2,

gdzie

-; = w(wH x)

jest rzutem prostopadłym wektora x ^{na wektor} w.

68 ^{A. K i e ł b a s i ń} s k i i K. Z i ^ę t a k

Unitarne przekształcenia Householdera nazywamy ^więc niekiedy przekształceniami odbić zwierciadlanych. Z tej interpretacji geometrycznej natychmiast wynika unitar-

ność przekształcenia P (tzn. Vx ^E K", llPxll ⁼ Will) oraz dalsze, ^łatwo sprawdzalne

własności:

pp= I, p = ^{p-1 =pH,}

_. I

x/ I I I I I -... I I

-... ... I

';-... . ..._

I

/ ^...

' I I I

I

I I

Rys. 1 (płaszczyzna rysunku ^rozpięta na wektorach w ^{i .X)}

Możemy również sprawdzić, że det(P) = - ^{1. Rozważmy} bowiem macierz

gdzie a ^{2 ,} a ^{3 , „.,} an ^są dowolnymi liniowo niezależnymi wektorami, rozpinającymi hiperpłaszczyznę prostopadłą do wektora w. Wówczas det(B) '# O.

Zbadajmy wyznacznik macierzy

PB= [Pw, Pa 2 , ••• ,Pan].

Ponieważ Pw = - w, Pai ⁼ ai, ^więc

det(PB) = det[-w, a ^{2 , •••} ,an]= -det(B).

Wobec tego, ^że det(PB) = det(P) · det(B), otrzymujemy det(P) = -1. Przekształ­

cenie Householdera może być oczywiście zapisane ^również za pomocą dowolnego wektora niezerowego u: ^{p = /-uTf~12 -UH, gdy w= u/llull ·c, lei = 1, c E K.}

2. Podstawowe zadanie (H) z przekształceniem Householdera. Rozważmy nastę­

pujące zadanie:

Dla danej pary niezerowych wektorów zespolonych a ⁱ e należących do K" chcemy wyznaczyć przekształcenie unitarne Householdera, które przeprowadza wektor a

na wektor Pa o kierunku wektora e, ^tzn.

(1) Pa = ek, lkl = llall/llell.

Ponieważ przekształcenie p jest hermitowskie, więc wyrażenie "QHp(i = aHek musi

być rzeczywiste. Okazuje się (zob. [6], str. 59), że jest to zarazem warunek dostateczny istnienia przekształcenia spełniającego (1.1) i (I). Zatem należy przyjąć

(2) k=

eHa 11-a11

I

± !(?Fial llell ' llall

± -llel! '

gdy

gdy Z wzoru (1) wynika

Niech (3) (4)

w= (a-ek)/lla-ekll.

u= a-ek,

R = I lull2 /2.

Wówczas przekształcenie (1.1) możemy zapisać w postaci

(5) P = 1-uRu .

^-+H

Parametr k określiliśmy z dokładnością do znaku. Znak ten wybieramy tak, by wektor u miał największą długość. Łatwo sprawdzić, że

(6)

Wobec tego we wzorach (2) wybieramy znak minus. A więc z (2), (3), (6) mamy

I-

^{-eu-a llal 1 a -+H-+}

U = a+e 1euaT ^{11e11 • g Y} e a# o,

-+ -+

llall d ~H-+ Q

a+e llell' ^{g Y} e a= ' (7)

(8) R = lla112+ leHalllallfllell.

Taki wybór znaku uzasadnimy później. Na razie ograniczymy się do uwagi, że

zapewnia on numeryczną poprawność (zob. [3]) najważniejszych algorytmów ko-

rzystających z przekształcenia Householdera.

z (2), (4), (6) otrzymujemy dla przypadku rzeczywistych wektorów a, e ^nastę­

pujące wzory na rzeczywiste przekształcenie Householdera:

(9) (10)

__ { a+esign(era)llall/llell, ^gdy

u - a+-e11a11111-e11, ^gdy

R = I lall2+1-era1 llall /Ile! 1.

era ^=I= O,

-era ⁼ o,

Konstrukcja rzeczywistego przekształcenia Householdera w rozważanym zadaniu ma prostą interpretację geometryczną. Mianowicie, wektor u ^leży na dwusiecznej

kąta między wektorami a ^{i -} e ^lub a ^{i e>} (patrz rys. 2).

70 ^{A. K ie ł b as i ń} s ki i K. Z i ^ę t a k

Dla powiększenia dokładności przekształcenia, wybieramy dwusieczną mniejszego

kąta (spośród kątów między wektorami a ⁱ e ^oraz a ⁱ -e). Wówczas kierunek od- powiedniej dwusiecznej jest „lepiej wyznaczony".

Podstawowym zadaniem korzystającym z przekształcenia Householdera jest wyznaczenie wektora c = Pb dla danego wektora bi danych wektorów a, e, ^defi-

/

,,.,,,.,,,.,,,.

,,.,,,.

_...,. -..„ ... __..

~,=a ·I e·k --- ----

Rys. 2 (płaszczyzna rysunku ^rozpięta na wektorach a ^{i e)}

niujących przekształcenie P. Zadanie to nazywać będziemy zadaniem (H). W dalszej

części pracy wykażemy, że wyznaczenie wektora c ^{może być} wykonane w sposób numerycznie poprawny, gdy przekształcenie P jest określone zależnościami (5), (7), (8). Chcąc wyznaczyć obraz c dowolnego wektora zespolonego b, nie musimy w sposób efektywny konstruować macierzy przekształcenia (5), ^gdyż

(11) c=Ph=b-u,RuHb. ~ _,. ( 1 _. _,)

Wygodnie jest zatem reprezentować przekształcenie P współrzędnymi wektora u

i liczbą R. Reprezentacja P zajmuje wtedy tylko ^około n miejsc (liczb rzeczywistych lub zespolonych) w ^pamięci maszyny (macierz P zajęłaby około n ² /2 miejsc) ^zaś koszt otrzymania c ze wzoru (11) jest ^rzędu tylko 2n ^mnożeń. (Zamiast n ² przy

posługiwaniu się pełną macierzą P~)

3. Zespolona arytmetyka zmiennopozycyjna cfl, cfl 2 , cfl. Zakładamy, że dyspo- nujemy dwójkową !-cyfrową arytmetyką zmiennopozycyjną f1 (zob. np. [5], str. 16,

[6], str. 110) ^określoną tak, ^że obliczone wyniki operacji arytmetycznych i pierwiast- kowania mogą być zapisane w postaci:

fl(xOy) = (xOy)(l-s), ^lei~ 2-t,

n(y':X) = Jlx (l - l5), I ól ^{~ p · 2-',}

gdzie O oznacza dowolny z operatorów ±, x, /,liczba p zaś należy na ogół do prze-

działu [I, 2]. Będziemy stosować relację ~ (wprowadzoną w [2] lub w [4]), co

. 1

pozwala na pomijanie w oszacowaniach błędów mniej istotnych składników rzędu

2- ^{21 ,} 2- ^{3 t,} itd. w obecności składników rzędu 2-'.

Niech zk = xk+ fyk (xb Yk rzeczywiste, f-jednostka urojona). Wówczas do- dawanie w arytmetyce cfl jest wykonane zgodnie z ^regułą

cfl(z1 +z2) ⁼ fl(x1 +x2)+ffl(Y1 +Y2) =

= (x1 +x2)0-e1)+i{y1 + Y2)(l -e2) = (z1 +z2)(l-s3), gdzie lsil ~ 2-t (s1, ^{s 2 -} rzeczywiste, s3 na ^ogół zespolona).

Podobnie wykonywane są w cfl pozostałe podstawowe operacje. Mianowicie,

łatwo pokazać, że

(I) cfl(z1/z2) = (z1/z 2 )(l-<p), cfl(lz1 ^{I) =} lz1 I ^(1-P),

cfl(rz1) ⁼ (rz1 ^{)(1- 6),}

lsl ~ (1 + J/2)2-t,

l'PI ^~ 1 (4 + y2)2-t,

IPI ^~ 1 (p+2,25)2-t'

161 ^~ 1 2-t,

gdzie r ^jest liczbą rzeczywistą, a lz 1 I, z1 /z 2 obliczamy według następujących wzorów (por. [7]):

I ^{0, gdy X1 = X2 = 0,}

lz1 ^{1 =} lx1IJ/I+(Y1/x1)2, gdy lx1I ~ IY1I, IY1IYl+(x1/Y1) ^{2 ,} gdy lx1I < IY1I,

Algorytmy te ^są „bezpieczne". ^Pozwalają na uniknięcie niedomiaru lub nadmiaru;

o ile dane i wyniki są dobrze reprezentowane w cfl. Dla przykładu podamy analizę numeryczną algorytmu dzielenia. Ze względu na „symetrię" algorytmu wystarczy

zbadać jeden przypadek. Np. lx 2 1 ^{~ IY} 2 1. ^Niech

Wówczas

(2) z= ^{cfl(~) =} X1+Y1h(l-01) (1- 03 )+t Y1 -Xi h(l- 04) (l- Os), Z2 X2 + Y2 h (1- 02) X2 + ^{Y2 h (I - 02)}

gdzie

1611, 1621, 1641~2.2-t, 1631, lósl ^{~ 3.2-t.}

1 1

Ponieważ sign(x 2) = sign(y 2 h), a ponadto ly 2 h/ml ^~ -!-, ^więc

x 2 + ^{y 2 h(I -} b,) = m( ^{1- y~h b,) = m(l-tp 2 ), l'1'2l ~ 2-'.}

72 ^A. Kiełbas iński i K. Z ^{i .ę} tak Zatem z (2) mamy

Z= _!_{[x1 +Y1h{l-'51)]{l-<p3)+f[y1-X1h(l-64)](l-<ps)}, l'Pil ~ 4.2-t.

m ^I

(3) Niech

Wówczas z (3) otrzymujemy

lz-wl ² lzl2

lz-wl ² lzl ²

Stąd, wykorzystując nierówność l.Z-zl ^~ lz-wl ⁺ lw-zl, otrzymujemy

co należało pokazać.

Zbadamy teraz błąd wytworzony przy obliczaniu w arytmetyce cfl iloczynu ska- larnego dwóch wektorów w ⁱ z ^o składowych wi, zi. Niech

(4)

Łatwo pokazać, że (por. (5], str. 30 i następne)

n n

(5) s ^{= cfl} (I: ^{ziwi) =} [I: ^z 1 ^{wj(I-~i)] (1- ~),}

j= 1 j= 1

l~jl ~ (n-j+ 1 + y2)2-i, ^{1~1 ~} 2-t.

1

A więc obliczony wynik sjest zaokrąglonym iloczynem skalarnym nieco zniekształ­

conych wektorów w ^lub z, ^gdyż czynnik (1- ^~ 1 ) możemy uważać za zaburzenie ^współ­

rzędnej w ₁ lub zi. Dla dowolnych wektorów w ⁱ z ^{nie ·} możemy twierdzić, że iloczyn skalarny jest obliczony z małym błędem względnym. Dużą dokładność względną obliczania wHz otrzymujemy w szczególnych przypadkach, np. gdy w ⁼ z. ^Miano-

wicie, (6)

n

cfl{z ⁸ z) ⁼ n(I: ^((rez 1) ² +(imzj)iJ) = (z ⁸ z)(I-<p),

J=l

l'PI ^~ (n+ 02-'.

I

Zajmiemy się teraz arytmetyką cfl 2 i cfl. Zakładamy, że czytelnik zna zasady arytmetyki fl 2 (zob. [5], str. 37) i fi (zob. [l]). Niech

n n

u= _i= L ₁ ^{ui =} L _{j= 1} ^(aj+ibi).

Wówczas w arytmetyce cf1 2 mamy

n n n n

(7) u = cfl2 ('L (aj+ ibi)) = fl2 (L ^{aj) +ifl2 (L>i) =} [L zj(1-1Pi)] (1-?p),

j=l j=l j=l i=l

11Pil ^~ Hn-j+ I)2- ² t, 111'1 ^~ 2-t.

1 1

Dla arytmetyki cfi otrzymujemy

n n

(8) u = cti(L ^{(ai+ibi)) =} l2:ZAI -?pi)] (1-?po),

j= 1 i= 1

Ostatecznie, łatwo pokazać, że dla iloczynu (4) obliczonego w arytmetyce cfl ₂ i cfl

otrzymujemy zależności (por. (5), (7), (8)):

n n

(9) s ^{= cf1 2} (L ^{z/wi) =} (L zi»HI-<pi)) (I-<p),

j=l 1

l<pil ~ t(n-j+I+ y2)2- ^{2 t,} !<pl~ 2-t,

1 1

n n

(10) s ⁼ cti(L ^{ziwi) =} [L ziwil-1Pi)] ^(1-?p),

j=I 1

Natomiast, ^jeśli z =w, to mamy (por. (6))

(11) cf1 ₂ (Z ⁸ z) = (Z ⁸ z)(l - a), lal ^~ 2-t,

1

(12) cfl(z ^{8 z) = (z 8 z)(l -?p), 11J'I ~ 3 • 2-t.}

1

Wobec tego kwadrat normy wektora z, ^a więc również jego norma, jest obliczany w arytmetykach cfl 2 , cfI z dużą dokładnością względną, niezależną od wymiaru n przestrzeni.

Pojęcie arytmetyki fl _{2 ,} a więc i cfl _{2 ,} nie jest ściśle określone. Najistotniejszą cechą

arytmetyk objętych tą nazwą jest opisana powyżej zdolność kumulowania iloczynów skalarnych na rejestrze podwójnej precyzji. W niektórych emc jest to jednak kosz- towne, ^więc stosuje się kumulację jedynie w przypadku sum (iloczynów skalarnych) wielkiej liczby składników. Jeśli dostęp do rejestru podwójnej precyzji nie jest kosz- towny, to możemy wykorzystywać go ^również do kumulacji nawet sum kilku ^skład­

ników, do operacji dzielenia z dzielną daną w podwójnej precyzji oraz do pierwiast-

kowania argumentu danego w podwójnej precyzji (por. [5], str. 38). ^Będziemy w ta-

kim przypadku mówili o ^pełnej arytmetyce ^ą ( ^cfl 2).

74 ^{A. K ie ł b as i ń s k i i K. Z i ę tak}

4. Zespolone zadanie (H) dla e ⁼ e1. ^{Niech aj, bj,} Uj oznaczają składowe wekto- rów a, b, u. Zakładamy, że składowe a ⁱ b nie są obarczone żadnym błędem i są równe swoim reprezentacjom w arytmetyce cfl. Niech

e = ei = (1, o, o, ... , o)r.

Jest to przypadek o szczególnie ^ważnych zastosowaniach. Wzory (2. 7), (2.8) ^można obecnie przekształcić do postaci (zob. [6], str. 308)

gdzie

R = A+S, u 1 = a 1 (I+ ~ ), u ₁ = a ₁ U= ^{2, 3, ... ,n),}

n

A = L _{i= I} ^{la1l 2 , s = v1a~l 2 A.}

Wyznaczając R i u1 z powyższych wzorów (zamiast z zależności

R = 11011 ² +Ja,111011, u, = a, (I ^{+ 1} 1 ^~; :)i

wykonujemy o jedno pierwiastkowanie mniej.

Postaramy ^{się podać} teraz realistyczne oszacowania ^błędów wytworzonych przy numerycznej realizacji zespolonego zadania (H). Przez F będziemy oznaczać nume- rycznie obliczoną wartość wyrażenia F (w pseudo-algolu możemy więc zapisać

F: ^{= F), a q;2-t} oznacza oszacowanie błędu cp (tzn. lcpl ^{~ q;2-t).}

I

Niech ex, {J, (}, µ oznaczają błędy względne, z jakimi są wyznaczone w danej arytmetyce wielkości A: ^{S, R, u 1 , tzn.}

(1) A-=A(ł-cx), S=S(l-cr), R=R(ł-e), il1=u1(1-µ),

v ⁼ [u1 , u2, ... , unl·

Wektor v jest identyczny z numerycznie wyznaczonym wektorem~. Zbadajmy związ­

ki zachodzące pomiędzy błędami cx, {J, ^(!, µ. Pozwoli nam to lepiej oszacować normę

macierzy ^błędów

(2) B = P-P, ^{P- = I _. -v-v} _R ^{1 -u}

niż gdybyśmy próbowali od razu ^korzystać ze związków pomiędzy oszacowaniami

a, a, e, µ.Łatwo sprawdzić, że

(3)

gdzie

a ^błędy Ll 1 , Ll 2 , Ll 3 mają oszacowania na poziomie 2- ^{2 t. Błędy} ex, e 1 , e 2 , e 4, ^u, e ^są

w istocie wyrażeniami rzeczywistymi, chociaż rozważamy tutaj zespolony przypadek zadania (H).

Zbadamy teraz, jak zaburzenia wektora u ^{i liczby R wpływają na} przekształcenie P. Z (ł) i (2) łatwo otrzymamy następujące wyrażenie dla B:

(4) B ₌ R -uu } ( ^--H e1+ue1u1µ+e1u U1µ-e1e1U1U1µµj, ^{-+-+T- - - -H - -T - ;;i} gdzie

ei = e+e2+e3+ ... = e+Ll4, ILl41 ^~ 2- 2 ^t,

u 1 , µ, ^są to liczby zespolone ^sprzężone z liczbani u _{1 ,} µ. Macierz B jest 1 więc hermitow- ska i ma następującą budowę:

Skorzystamy teraz z tej szczególnej postaci B. Zauważmy, że pomijając elementy pierwszego wiersza, wszystkie kolumny tej macierzy są równoległe do wektora (O, u2, u3, ... , un)T. Jeśli więc pomnożymy lewostronnie macierz B przez macierz odbicia zwierciadlanego G, przeprowadzającego (n- !)-wymiarowy wektor [u 2 , u 3 ,

... ' unY na wektor o kierunku wektora e2' to otrzymamy macierz GB, ^mającą nie- zmieniony pierwszy wiersz, ^zaś wiersz drugi postaci

! (Uu1(µ-e1) - UU2e1 ... - UUnei), gdzie

n

u= (L ₎₌₂ ^{1uj1 2} f' ^{2 •}

Następne wiersze zawierają wyłącznie elementy zerowe. Zatem wszystkie wiersze macierzy GB, z pominięciem elementów pierwszej kolumny, są bądź równoległe do wektora (O, u ^{2 ,} u 3 • ••• , unY, bądź składają się z samych zer. Zatem macierz GBGH ma wszystkie elementy równe zeru z ^wyjątkiem czterech elementów w lewym górnym rogu. Ten niezerowy blok ma ^postać

(6) - 1 (µ+P,-µP,-e1)U1U1 (µ-e1)U1U)

C-R (µ-e1)ii1U -e1U2 .

76 ^{A. K i e ł b a s i ń s} k i i K. Z i ^ę t a k

Oczywiście norma spektralna macierzy B jest równa normie spektralnej macierzy C.

Aby wyznaczyć tę normę, musimy znaleźć największy moduł wartości własnej ma- cierzy C. Wielomian charakterystyczny ma ^postać

(7) A ² + A.[2e1 ^-(I +h)(µ+µ-µµ)]-(l-h ² )(l -e1)µµ = o.

Widzimy ^więc (zob. (3)), że maksymalny moduł wartości własnej jest równy

IAlmax = ł{l2e1 -(1 +h)(µ+/,t-µ/,t)I+

+ y[2e1 ^-(1 +h)(µ+µ-µµ)]2+4(1-h ² )(l-e1)µµ J.

Skorzystamy teraz ze ^związków (3)

(8) IAl-lłcx+e1h - l+h ^+e1 - e2+e4 ^-(1 + h)e3+e3+L1sl 2 +

-. /(tcx+e1 +h (l h) e3+e3+L1s)2

+V l+h +e1-e2+e4- + 2 +

gdzie L1 i są rzędu wielkości 2- ² r.

Oznaczamy oszacowanie prawej strony powyższej równości przez y2-t. ^Zatem

dla & (mnożnik z oszacowania ^błędu ex) i y otrzymujemy dla p = 1 (zob. (3.6), (?.11), (3.12), (3)), ^wartości przedstawione w ^poniższej tabelce:

~

ot

'Y

cfl n+l

łn+9+ v<tn+9)2+(łn+6) ² ~

~ (1 + y'2) ^(łn+9)

Z powyższych rozważań mamy oszacowanie postaci

(6) llP-Plb ⁼ llBlb ^~ 2-'y.

cii ^cf}z

3

24 23,5

Niech g oznacza numerycznie wyznaczony obraz wektora b. Porównajmy go z wektorem 7: = Pb:

(:10) llg-cll ^~ llg-Phll + llPb-Pbll ^~ lig-Phil +1-rylfbll.

Oszacujemy teraz pierwszy ^składnik. Wektor g wyznaczamy w dwóch krokach

(~ob. (2.11))

-+Hb-+

- · - V

w.----, _R g ^:= h-vw.

Wobec tego z (3.1), (3.5), (3.9), (3.10) otrzymujemy

C.· ~

(11) vH(f-D1)b

w = - - - -

R

gdzie

D 1 = diag(di), D 2 = diag(Bi), D 3 = diag(ai),

a błędy (jb {Ji, rx.i, mają oszacowania z mnożnikami ()~, ~, &i odpowiednio równymi:

cfl cti cfh

n-j+3+ y.f 4+y.f 2

&i = 1 + Jf2, ^{Pi = 1 dla każdej z} rozważanych arytmetyk. Zatem stąd oraz z (2), (3) mamy następujące oszacowanie:

(12)

~ (2Ó1 +3fl1 +2&1)11hll2-t = i11b'112-t.

1

Ostatecznie więc z (10), (12) wynika, że

llg-cll ~ ₁ (y+ d)llbll2-' = Kwllc112-t,

przy czym (przyjęliśmy p = 1) cfl 2n+ 11+3 y'.f

3,2n+36

cii 13+4y.f

42,6

9+2V.f 35

Widzimy więc, że sądząc z oszacowań błędu, arytmetyka cfI jest prawie równie dobra, jak arytmetyka cf1 _{2 •}

W obu arytmetykach wektor c jest wyznaczany z dużą dokładnością względną, niezależną od wymiaru n.

5. Rzeczywiste zadanie H dla e ⁼ ep ^Jeśli e = e ^{1 ,} to wzory (2.9), (2.10) przyj-

mują postać

gdzie

gdy ll1 ~o,

gdy a 1 <O, A = lla112, N = llall .

ui = ai (j = 2, 3, ... , n),

J. H. Wilkinson przedstawia w [6], str. 148 i następne, analizę błędu wytworzonego przy numerycznej realizacji w pełnej arytmetyce fl 2 rzeczywistego zadania H dla

e ⁼ e ^{1 •} w niniejszej pracy podamy oszacowania potencjalnie nieco lepsze od osza- cowania otrzymanego przez Wilkinsona (zob. ^również [l]). Nadal stosujemy ozna- czenia wprowadzone w § 4. Niech

A= ^A(l-a), N= N(l-v), u 1 = u 1 (1-µ), R ^{= R(l-e).}

78 ^{A. K i e ł b a s i ń s k i i K. Z i ę t a k}

Wprowadzimy obecnie zależności pomiędzy błędami ac, „, ^{e, µ.} Łatwo sprawdzić, że

1-'P = ^~ _(1-e1), ^{le1I ~ p2-•,}

1-µ=(1- _a +·„N( )N)(t-e2), le21~2-•,

1 s1gn a1

1- - e =(i- ocA+(v+s3-e3v)la1IN)(l- ) ^{A+la1IN E4'} I ^{E3, E4} I I I~ ^~ ₂ ^{-. i}

Stąd wynika, ^że

„ ^{= łcx+e1 +L'.11,}

(1) µ = ^łoc+e1 1 ^{+h +e2} +L'.12,

n = a+ ^{(łoc+e1 +e3)h L1}

t:' 1 +h +e4+ 3,

gdzie

la1I

h = lllill '

a błędy L1i mają oszacowania na poziomie 2-2t. Chcemy podkreślić, że rozważamy tu obliczenie R w „niepełnej" arytmetyce f1 2 , w następujący sposób

(2) R ⁼ fl[fl(fl2(llal12))+ yfl{fl2(llal12)) la1I].

Pokażemy teraz, jak zniekształcenie wektora u i parametru ii. ^{wpływa na prze-}

kształcenie P. Niech tak jak w § 4

llBll ⁼ llP-Pll ^~ 2-'y.

Podobnie jak w poprzednim punkcie

B ₌ R 1 ( ... ... ^-uu T ei +e1U U1µ+ue1 U1µ-e1e1U1µ ' ... ...T .... ... T

T 2 2) gdzie

ei ⁼ e+e ² +e ³ + ^{„. =} e+L'.14, IL'.141 ^~ 2- ^2'.

1

Błąd a dla arytmetyki fl 2 i fl ma takie samo oszacowanie jak w § 4 dla arytmetyki

cfl 2 i cfl, a dla arytmetyki fl możemy pokazać, że lal ^~ n2-'. Odpowiednia tabelka dla y wygląda więc następująco: 1

fi fi

y

łn+4+ Jl<tn+4) ² + (łn+2) ² < (1 + y2)(łn+4) 12 9,7

Uwzględniając analogiczne zmiany w oszacowaniach elementów macierzy Di otrzy- mujemy dla błędu otrzymanego obrazu g wektora b oszacowanie:

llg-cll ^~ (y+ó)llbi12-' = Kwllcll2-',

1

przy czym

fl fi

2n+7 13 9

3,2n+17 25 18,7

Otrzymane wartości wskaźnika kumulacji algorytmu, Kw, ^są dla arytmetyk fi ^{i fl} 2 większe od ^wartości podanych w [l] (porównaj ^również [6], str. 154). Jak ^już wspomi-

naliśmy, jest to ^związane m.in. ze sposobem obliczania ^wartości parametru R (zob. (2)). Ponadto, ze ^względu na ograniczenie ^zasięgu stosowania podwójnej precyzji w artmetyce fl 2 jedynie do kumulowania sum wielu składników, otrzymu- jemy większe oszacowania błędu wytworzonego przy obliczaniu Pb. Mianowicie, pokazaliśmy, że mnożnik J ⁼ 9, gdy tymczasem Wilkinson, zakładając pełne wy- korzystanie podwójnej precyzji w arytmetyce fl 2 , otrzymuje 3 ⁼ 3,35. Z tych samych

względów istnieje możliwość obniżania oszacowań otrzymanych dla zespolonego zadania H, ^gdybyśmy rozszerzyli zakres stosowania podwójnej precyzji w arytme- tyce cfl _{2 •}

Tak więc pokazaliśmy, że zespolone i rzeczywiste zadanie (H) dla e ⁼ e ^{1 może}

być realizowane w sposób numerycznie poprawny w każdej z rozważanych arytme- tyk.

6. Numeryczna reaJizacja ogólnego zadania (H). Niech Pa oznacza przekształcenie

Householdera przeprowadzające wektor a na wektor o kierunku wektora e. ^Wówczas

Pa jest reprezentowane przez wektor u(a) i parametr Ra (zob. (2.5)-(2.7)).

Wprowadzimy pomocnicze oznaczenia: v- numerycznie wyznaczony wektor u(a), 1- zaburzony Wektor ll, taki, by leHll)obl = eH/: P:- numerycznie wyzna- czane przekształcenie Householdera reprezentowane przez v ^{i Ra, g - numerycz-}

nie wyznaczony wektor Pab.

W przeprowadzonej przez nas analizie numerycznej ogólnego zadania (szczegóły

pomijamy) otrzymaliśmy następujące oszacowania

111-all ^~ Kd 11a112-t, llv-u(f)ll ^~ µIlu (f)ll 2-t,

1 1

llPa-P1ll ^~ ,Y2-t, llg-P1bll ~ Kw1lbil2-t = Kwllcll2-t,

1 1 1

gdzie c ^{= P);, P} 1 oznacza przekształcenie Householdera skonstruowane dla wektora

fi e, ^a Kd, Kw, y, P, są odpowiednio równe:

fl fi f h ^{cfl cti cfl2}

" ^7n+2,8 4,9 3,5 2,ln+ 11,6 15,8 9,6

µ

A

6,8n+ 13 33,4 22,8 14,5n+56 91,6 50,2

y

Kd o o o n+y2 2+y2 ^<2-f)

Kw ^{8,8n+20 47 32} I 16,5n+71 110 62

80 ^{A. K i e ł b a s i ń} s k i i K. Z i ^ę t a k

Widzimy więc, że w przypadku zespolonym otrzymany obraz g wektora b ^dla

transformacji Pa może być interpretowany jako trochę zaburzony obraz wektora b

dla transformacji P1, przy czym wektor f- a jest małym zaburzeniem wektora a.

W tym sensie realizację ogólnego zadania (H) możemy uważać za numerycznie poprawną. w przypadku rzeczywistym f ⁼ a, w przypadku zespolonym zaś dla arytmetyki cfl 2 zaburzenia/_:_ a ^są na poziomie 2- ² 'llal 1. W trakcie analizy numerycznej zespolonego zadania (H) okazało się, że o dokładności obliczonego wektora P

b

decyduje przede wszystkim dokładnoŚĆ, Z jaką oblicza ^się iloczyn Skalarny eHa.

o wrażliwości wektora u(a) na zaburzenia wektora a decyduje bowiem wrażliwość

ilorazu eHa/leHal ^{na małe zmiany} składowych Wektora Q (a niekiedy ^{może być} Ona

duża).

Gdyby z ^jakichś powodów potrzebna była (również dla przypadku zespolonego) wysoka dokładność względna wektora g, należałoby wtedy ^stosować jedynie arytme-

tykę cfl _{2 •}

Bibliografia

[1] A. K ie ^ł bas i ^ń s k i, Algorytm sumowania z poprawkami i niektóre jego zastosowallia, Matematyka Stosowana 1 (1973), str. 23-41.

[2] - Allaliza numerycz!la algorytmu ortogonalizacji Grama-Schmidta, ibid. 2 (1974), str. 15-35.

[3] - Podstawowe ^pojęcia analizy ^błędów w metodach numerycznych algebry liniowej, ibid. 4 (1975), str. 5-27.

[4] J. M aj c h r o w s ka i A. S m o k t u n o w i cz, Analiza numeryczna algorytmu Ortegi- Householdera w dziedzinie zespolonej, ten tom, str. 55-66.

[5] J. H. ^Wił ki n son, Błędy zaokrągleń w procesach algebraicznych, Warszawa 1967.

[6] - A1ue6pau11.ec1<aR npo611eMa co6cm6e111,blx JHa11.euui'i, MocKBa 1970.

(7) J. H. W i I ki n son and C. Re i n s c h, Handbook for automatic computation, Springer-

Verlag, 1971, str. 359-371 (lub Numer. Math. 12 (1968), str. 349-368).