Procedura ortogonalizacji układu wektorów

ROCZNIKI POLSKIEGO TOWARZYSTWA MATEMATYCZNEGO Seria III: MATEMATYKA STOSOWANA VI (1976)

ANNA GÓRAJ i ANDRZEJ KIEŁBASIŃSKI (Warszawa)

Procedura ortogonalizacji ^układu wektorów

metodą Grama-Schmidta

1. ^Wstęp. Przedstawiamy procedury GSPO i GSPO 1, realizujące algorytm Grama-Schmidta z poprawianiem (por. [l]). Każda z nich ortonormalizuje układ r niezerowych wektorów, danych jako kolumny macierzy X (n x r ), n ;::;: r > O. W tym celu macierz X jest rozkładana na iloczyn macierzy Y (n x r) o zortonormalizowa- nych kolumnach i macierzy górnej trójkątnej C (r x r) o dodatnich elementach na ·

przekątnej. Algorytm GSPO realizuje obliczenia w standardowej arytmetyce fl (por. [3]), zaś GSPO 1 w fl (por. [2]). Testowanie przeprowadzono w Zakładzie

Obliczeń Numerycznych U.W. na maszynie GIER.

2. ^Pa~ametry procedur

r - wartość całkowita, liczba wektorów ortonormalizowanych,

n - wartość całkowita, wymiar przestrzeni euklidesowej wektorów, n ;::;: r, ro - wartość rzeczywista, najmniejsza liczba taka, ^że reprezentacja w fl liczby

1 +ro jest większa niż I (zazwyczaj ro = 2-t, t liczba cyfr mantysy binar- nej w fl),

X, Y - tablice rzeczywiste dwuwskaźnikowe o zakresie [l :n, I :r]. Przy ^wejściu X zawiera ^układ ortonormalizowanych wektorów. Przy ^wyjściu Y zawiera

układ ortonormalnych wektorów.

C- tablica rzeczywista dwuwskaźnikowa o zakresie co najmmeJ [l :r, 1 :r].

Przy wyjściu zawiera w górnym prawym trójkącie macierz trójkątną górną rozkładu X= Y · C z dodatnimi elementami na przekątnej.

Uwag a. Przy ^wywołaniu procedur możemy przyjąć parametr aktualny Y lub C - identyczny z parametrem aktualnym X. ^Wywołanie procedur niszczy wówczas

zawartość pierwotną X. Parametry aktualne Y i C muszą być oczywiście różne.

3. Metoda testowania. Procedury testowano ^głównie na macierzach X liczb pseudolosowych o rozkładzie jednostajnym z przedziału [ - 10, 10], porównując

z działaniem procedury GSP02, wykorzystującej kumulowanie iloczynów skalar- nych na 39-u bitach (w GIER fl liczba bitów mantysy t = 29). Stosując taką samą kumulację, obliczano residua względne (w normie euklidesowej):

W1 = llX-YCll/llXll,

Ogółem przebadano ^około 60 przypadków dla różnych n i r. Grube oszacowania teoretyczne {por. [1]) dla W

W

podaje tabelka I (dla arytmetyki GIER ^należy ·

położyć ro = 4

9):

Tabelka 1 fi

2(r+ yr). ^ro _{4 •} V~•

^{5 ·}

^I

1_ __

w_2 _~ _{____} n_·_v_r_·_ro ___ ---'---~~ y ^{r ·} r~--·--5-=-· '~--- _ _\

Uzyskane w eksperymencie ^wielkości WP okazały się od kilku do kilkudziesięciu

razy mniejsze od oszacowań teoretycznych. Większa dokładność arytmetyki fi

(GSPOl) w stosunku do fi (GSPO) dała się zauważyć dla dużych wartości n. Przy-

kładowo podajemy kilka maksymalnych wartości, osiąganych przez WP w ekspe- rymencie:

Tabelka 2

-~z;~;~~ów J ___ n __ !_J ^p GSPO ·- _ __ G~-~0_1 _J_ ^{Gs:o~- _}

2 35 ³⁵ I - ·- -1---1---1--- ^{J 1710-9}

2

5 400 4

2 ' - - - - ---

4. Szczegóły organizacyjne i zakres stosowalności procedur

4.1. Przy wywołaniu tablica X nie może zawierać wektorów zerowych, ^gdyż dla zabezpieczenia przed przekroczeniem zakresu liczb dobrze reprezentowalnych w arytmetyce fi (por. [3]) w procedurach wykonywane jest skalowanie kolumn do poziomu ^jedności. Nie ma natomiast innych ograniczeń danych. Tablica X ^może

zawierać nawet ^układ wektorów liniowo zależnych. Jakość otrzymanego ^rozkładu X = Y · C charakteryzuje a priori tabelka 1.

Uwag a. Liniowa ^zależność (w sensie arytmetyki fi) wektora x z k-tej kolumny X od „wcześniejszych" kolumn macierzy Yjest sygnalizowana ^relacją:

C[k, k] <ro· ll:Xll · ^g,

g = ^2k dla GSPO, g = 3yk dla GSPOl.

4.2. Procedury wykorzystują około r + 15 zmiennych roboczych. Realizacja pro-

cedur wymaga ponadto co ^najwyżej r(2 · n+r) (lub mniej) miejsc na tablice X, Y, C

w ^pamięci emc. Wywołanie procedur instrukcjami postaci: OSPO (r, n, ro, A, B, C),

lub OSPO (r, n, ro, B, B, C), lub GSPO (r, n, ro, A, B, A) pozwala zachować za-

Procedura ortogonalizacji

wektorów

Grama-Schmidta 19

wartość tablicy A nienaruszoną lub „zapisaną" czynnikiem ortonormalnym, lub czynnikiem trójkątnym.

4.3. Czas potrzebny do realizacji procedur nie może być ściśle określony z góry,

gdyż zależy od ilości „poprawień". Orientacyjna ^formuła asymptotyczna (dla do- statecznie dużych n i r): n · r

(M + ^S + ^{2 · T} + ^F) podaje czas minimalny (gdy nie ma żadnych poprawień) dla GSPO. (M - czas mnożenia, S - czas sumowania, T- pobrania elementu tablicy dwuwskaźnikowej, F- czas organizacji ^pętli for.) Dla GSPOl należy S zastąpić przez 3 · S. Przy poprawianiach czas może wzros-

nąć co najwyżej trzykrotnie.

5. Algorytmy

procedure GSPO (r, n, ro, X, Y, C);

value r, n, ro;

integer r, n ; real ro;

array X, Y, C;

begin

integer i,p, k, kl ,j;

real NI, N2, N3, eps, g, z;

array f[l :r-1];

ro : = 2 x ro x ro;

for k : = 1 step 1 until r do begin

kI:=k-1;

N2:= g:= O;

for i : = 1 step 1 until n do begin

z:= abs(X[i, k]);

if z > g then g : = z end;

g:= 2 tentier (ln(g)/0.6931);

for i : = I step 1 until n do begin

z:= ^{Y[i, k]} := X[i, k]/g;

N2:= N2+zxz end;

for p: = 1 step I until kl do C[p, k] :=O;

eps : = k x k x N2 x ro;

for j:= 1, j+ 1 while N2- NI > N2/4/k do begin

if j > I then

I

begin ' N2:= Nl;

if ^NI~ eps then Y[k, k] := sqrt(eps) end;

for p : = I step I until kl do begin

NI:= O;

for i : = 1 step 1 until n do NI:= Nl+Y[i,p]x.Y[i, k];

f[p] := Nl;

C[p, k]: = C[p, k]+ NI x g end;

Ni:= O;

for i : = I step I until n do begin

end end j;

N3:= O;

for p: = I step 1 until kl do N3 := N3+ Y[i,p] xf[p];

N3:= Y[i,k]:= Y[i,k]-N3;

NI:= Nl+N3xN3

NI:= sqrt(Nl);

C[k, k] :::;= NI xg;

for i: = I step 1 until n do Y[i, k]: = ^{Y[i, k]/ NI}

end k end GSPO;

procedure GSPOI (r, n, ro, X, Y, C);

value r, n, ro;

integer r, n ; real ro;

array X, Y, C;

begin

integer i,p, k, kl ,j;

real NI, N2, N3, eps, g, t, tl, z, v;

array /[l :r- I];

ro : = 5 x ro x ro;

for k : = 1 step I until r do begin

kl:=k-I;

N2:= g:= O;

\

Procedura ortogonalizacji ^układu wektorów ^metodą Grama-Schmidta 21 for i : = 1 step I until n do

begin

z:= abs(X[i, k]);

if z > g then g : = z end;

g: = 2 tentfer(ln(g)/0.6931);

for i : = 1 step I until n do begin

z:= Y[i, k] := X[i, k]/g;.

N2:= N2+zxz end;

for p : = I step I until k I do C[p, k] :=O;

eps : = k x N2 x ro;

for j: = 1, .i+ 1 while N2-Nl > N2/4/k do begin

if j > I then begin

N2:= Nl;

if Nl ~ eps then Y[k, k]: = sqrt(eps) end;

for p: = l step I until kl do begin

NI:= t:= O;

for i : = 1 step 1 until n do begin

V:=

Y[i, p]

Y[i, k];

z:= Nl+v;

t:= Nl-z+v+t;

Ni:= z end;

f[p] : = N 1 : = N 1 + ^t;

C[p,k]:= C[p,k]+Nlxg end;

N3:=tI:=O;

for i : = I step I until n do begin

Nl:=t:=O;

for p : = 1 step l until k I do begin

V :

= Y[i' p]

f[p];

z:= Nl+v;

t:= Nl-z+v+t;

NI:= z

'

end;

V:= ^{Y[i, k]:} = Y[i, k]- (NI +t);

v := vxv;

z:= NJ+v;

tl := N3-z+v+tl;

N3:= z end;

NI:= N3+tI end j;

Nl := sqrt(Nl);

C[k, k] := Nl xg;

for i : = I step I until n do Y[i, k]: = Y[i, k]/ NI end k

end GSPOI;

Cytowana literatura

[1] A.

Analiza numeryczna algorytmu Grama-Schmidta, Mat. Stos. 2 (1973), str. 15-35.

[2] . - Algorytm sumowania z poprawianiem ⁱ niektóre jego zastosowania, Mat. Stos. 1 (1973), str. 24-41.

[3] J. H. W i I ki n son,

w procesach algebraicznych, Warszawa 1967.