Za realne brojeve definiramo funkciju “apsolutna vrijednost”: Prelaskom na skup kompleksnih brojeva, javila se potreba za prirodnim poopćenjem te funkcije: Poanta takvih poopćenja je primjena funkcije na širem skupu ulaznih parametara, uz zadržavanje ponašanja na prvotnom (užem) skupu.
U ovom članku bavimo se matričnim poopćenjima uobičajenih funkcija koje smo navikli primjenjivati na brojevima. Dakle, definirat ćemo funkcije koje (neke) matrice preslikavaju u druge matrice, pri čemu se zadržavaju dobra svojstva odgovarajućih funkcija na skupu kompleksnih brojeva. Kao osnova članka korišten je diplomski rad N.Strabić [6], koji je baziran pretežno na knjizi N.Highama [1].

Za početak, potrebno je naglasiti kako matrične funkcije imaju širok spektar primjena u znanosti i inženjerstvu. Naprimjer, uz polarnu dekompoziciju matrice A vežemo matrični korijen, jer znamo da je P = \sqrt{A^{*}A}. Zanimljivo, i za računanje polarne dekompozicije i za računanje matričnog korijena korisno je definirati matričnu funkciju predznaka kojom se bavimo u nastavku ovog rada.

S funkcijom predznaka obično se susrećemo u ranom obrazovanju, jer je riječ o jednostavnoj, ali ipak “neobičnoj” funkciji (u smislu da nije definirana “običnom” formulom nego rastavljena na slučajeve): Primijetimo da vrijedi: Kao što smo vidjeli u (1), funkciju |\cdot| definiramo i za kompleksne brojeve, pa (2) možemo upotrijebiti kao definiciju funkcije \mathop{\text{sign}} x za sve kompleksne brojeve x različite od nule. No, takva funkcija predznaka neće više nužno biti realna, nego kompleksna, a \mathop{\text{sign}} x predstavlja točku na jediničnoj kružnici koja je najbliža zadanom kompleksnom broju x.

Budući da obično želimo zadržati realnost funkcije predznaka i za kompleksne argumente, funkciju definiramo u ovisnosti o realnom dijelu kompleksnog broja: pri čemu s \mathbb{I} označavamo skup imaginarnih brojeva, a s \Re x realni dio kompleksnog broja x. Možemo li sličnu generalizaciju napraviti i za matrice?

Nažalost, skok s brojeva na matrice bitno je složeniji nego skok s realnih na kompleksne brojeve. Naprimjer, množenje matrica bitno se razlikuje od množenja brojeva: za matrice A i B lako se može dogoditi da je definiran umnožak AB, ali ne i BA. Čak i kad su oba produkta definirana, oni ne moraju biti jednaki.

Kako bismo izbjegli zavrzlame oko toga kada je produkt dviju matrica definiran, u nastavku ćemo razmatrati kvadratne matrice istog reda.




Za proizvoljnu matricu A\in \mathbb{C}^{n \times n} i proizvoljni polinom p u varijabli t, matrica p(A) definira se prirodno, zamjenom t s A, što možemo izvesti, jer su za kvadratne matrice prirodno definirane nenegativne cijele potencije: pri čemu s I_{n} označavamo jediničnu matricu reda n. Ako je red matrice jasan iz konteksta, umjesto I_{n} pisat ćemo I.

U ovom poglavlju bavimo se općenitijim slučajem: za kompleksnu funkciju f i matricu A \in \mathbb{C}^{n \times n} definiramo matricu f(A). Prvo ćemo dati dvije moguće definicije matričnih funkcija, a zatim i iskaz teorema koji upućuje na istovjetnost tih dviju definicija.


Prva definicija matričnih funkcija koju ćemo razmotriti bazira se na Jordanovoj kanonskoj formi matrice. Dat ćemo kraći pregled potrebnih pojmova i rezultata, ali korištene rezultate nećemo dokazivati. Zainteresirani čitatelji dokaze mogu naći u mnogim knjigama koje se bave osnovama teorije matrica, npr.u [3], [4], [5] i [2].


Za matricu A\in \mathbb{C}^{n \times n} kažemo da je singularna ako je \det A = 0, što je ekvivalentno tome da matrica nema inverz, tj.da ne postoji matrica B\in \mathbb{C}^{n \times n} takva da je AB = I (ili, ekvivalentno, BA = I). Ako matrica A ima inverz (kažemo da je regularna ili invertibilna), taj inverz je jedinstven i označavamo ga s A^{-1}.

Skalar \lambda \in \mathbb{C} zovemo svojstvena vrijednost matrice A\in \mathbb{C}^{n \times n} ako postoji vektor v \in \mathbb{C}^{n} \setminus \lbrace 0\rbrace takav da je Vektor v se zove svojstveni vektor matrice A pridružen svojstvenoj vrijednosti \lambda. Skup svih svojstvenih vrijednosti matrice A nazivamo njenim spektrom, u oznaci \sigma(A). Vrijedi da je \lambda \in \sigma(A) ako i samo ako je A - \lambda I singularna matrica, odnosno ako i samo ako je \det (A - \lambda I) = 0.

Polinom k_{A}(t) := \det(A - tI) nazivamo karakteristični polinom matrice A i obično ga označavamo samo s k(t). Očito, k(t) je polinom n-tog stupnja, a njegove nultočke su upravo svojstvene vrijednosti matrice. Vrijedi i k(A) = 0.

Za matricu B\in \mathbb{C}^{n \times n} kažemo da je slična matrici A\in \mathbb{C}^{n \times n} ako postoji (regularna) matrica S\in \mathbb{C}^{n \times n} takva da je Lako se pokaže da je matrična sličnost relacija ekvivalencije.
Jordanov blok J_{k}(\lambda) je gornjetrokutasta1 matrica reda k oblika Jordanova matrica J\in \mathbb{C}^{n \times n} je blok-dijagonalna matrica oblika

J = \mathop{\text{diag}}(J_{k_{1}}(\lambda_{1}),\dots, J_{k_{s}}(\lambda_{s})) = \begin{bmatrix} J_{k_{1}}(\lambda_{1}) \\ && \ddots \\ &&& J_{k_{s}}(\lambda_{s}) \end{bmatrix},

gdje je k_{1} + k_{2} + \dots + k_{s} = n red matrice. Redovi k_{i} Jordanovih blokova ne moraju biti međusobno različiti, kao niti njihove svojstvene vrijednosti \lambda_{i}.

Teorem 1. (Jordanova kanonska forma) Za svaku matricu A\in \mathbb{C}^{n \times n} postoji Jordanova kanonska forma J_{A}, odnosno matrica oblika

J_{A} := \mathop{\text{diag}}(J_{k_{1}}(\lambda_{1}), J_{k_{2}}(\lambda_{2}),\dots, J_{k_{s}}(\lambda_{s}))

slična matrici A, pri čemu su \lambda_{i} sve svojstvene vrijednosti matrice A.

Matrica J_{A} je jedinstvena, do na poredak blokova (koji možemo proizvoljno izabrati odgovarajućom promjenom u matrici sličnosti). Matrica sličnosti nije jedinstvena2.

Red najvećeg Jordanova bloka koji odgovara svojstvenoj vrijednosti \lambda, dakle, broj \max\lbrace k_{i}:\ \lambda_{i} = \lambda\rbrace, zovemo indeks svojstvene vrijednosti \lambda. Broj blokova koji odgovaraju određenoj svojstvenoj vrijednosti naziva se geometrijska kratnost. Zbroj redova svih tih blokova, tj.broj pojavljivanja broja \lambda na dijagonali Jordanove forme, zovemo algebarska kratnost te svojstvene vrijednosti.

Jordanova kanonska forma je vrlo korisna za teorijska razmatranja, jer mnoga prirodna matrična svojstva i operacije čini preglednima. Naprimjer, iz same forme lako vidimo rang matrice (broj elemenata dijagonale različitih od nule), svojstvene vrijednosti i njihove pripadne kratnosti. Potenciranje matrice svodi se na potenciranje matrica J_{k_{i}}(\lambda_{i}), za što postoji jednostavna egzaktna formula, računanje determinante svodi se na množenje dijagonalnih elemenata matrice J_{A} itd.

U praksi se, nažalost, Jordanova forma mora vrlo oprezno koristiti zbog numeričke nestabilnosti. Male greške, nastale zaokruživanjem zbog konačnosti računalne aritmetike ili nepreciznosti mjerenja, mogu dovesti do bitne promjene tipa Jordanove forme — raspad velikog Jordanova bloka u niz manjih, raspad višestruke svojstvene vrijednosti u nekoliko bliskih jednostrukih itd.

2.1.2Definicija matrične funkcije

Označimo sve različite svojstvene vrijednosti matrice A s \lambda_{1}, \dots, \lambda_{s} i njihove pripadne indekse s n_{1},\dots, n_{s}.

Za kompleksnu funkciju f kažemo da je definirana na spektru od A ako su definirane sve vrijednosti

f^{(j)}(\lambda_{i}), \quad j=0,\dots,n_{i}-1, \quad i=1,\dots,s.

Primijetimo da za definiranost na spektru nije dovoljno definirati samo vrijednosti funkcije za svojstvene vrijednosti, nego i vrijednosti odgovarajućeg broja derivacija funkcije f u svojstvenim vrijednostima matrice.

Definicija 2. [Matrična funkcija preko Jordanove forme] Neka je f kompleksna funkcija definirana na spektru matrice A\in \mathbb{C}^{n \times n}, J Jordanova kanonska forma matrice A i S regularna matrica takva da je A = S^{-1}JS. Tada je

(4)

f(A) := S^{-1}f(J)S = S^{-1} \mathop{\text{diag}}(f(J_{k_{1}}(\lambda_{1})),\dots,f(J_{k_{s}}(\lambda_{s}))) S,

pri čemu je

(5)

f(J_{k_{i}}(\lambda_{i})) := \begin{bmatrix} f(\lambda_{i})& \displaystyle \frac{f'(\lambda_{i})}{1}& \displaystyle \frac{f''(\lambda_{i})}{2!}& \displaystyle \dots& \displaystyle \frac{f^{(k_{i}-1)}(\lambda_{i})}{(k_{i}-1)!} \\ & \displaystyle f(\lambda_{i})& \displaystyle \frac{f'(\lambda_{i})}{1}& \displaystyle \dots& \displaystyle \frac{f^{(k_{i}-2)}(\lambda_{i})}{(k_{i}-2)!} \\ && \displaystyle f(\lambda_{i})& \displaystyle \dots& \displaystyle \frac{f^{(k_{i}-3)}(\lambda_{i})}{(k_{i}-3)!} \\ &&& \displaystyle \ddots & \displaystyle \vdots \\ &&&& \displaystyle f(\lambda_{i}) \end{bmatrix}.

Prisjetimo se, matrica J jedinstvena je do na poredak dijagonalnih blokova, ali matrica sličnosti S nije jedinstvena. Zbog toga, potrebno je dokazati ispravnost prethodne definicije, tj.pokazati da ona ne ovisi o izboru matrice S.

Propozicija 3. Definicija 2 ne ovisi o izboru Jordanove kanonske forme matrice A.

Skica dokaza propozicije dostupna je u Problemu 1.1 u [1].

Primjer 4. Neka su zadane matrica

A = \begin{bmatrix} 2 & 3 & 1 \\ 0 & 0 & 0 \\ -4 & -6 & -2 \end{bmatrix}

i funkcija

f(t) = \sin t - e^{t}.

Želimo izračunati f(A).

Najprije tražimo Jordanovu kanonsku formu, tj.matrice J i S takve da je A = S^{-1}JS. Jednostavno se vidi da je

S^{-1} = \begin{bmatrix} -1 & -1/2 & -3 \\ 0 & 0 & 2 \\ 2 & 0 & 0 \end{bmatrix},\ S = \begin{bmatrix} 0 & 0 & 1/2 \\ -2 & -3 & -1 \\ 0 & 1/2 & 0 \end{bmatrix} \text{ i } J = \begin{bmatrix} 0 & 1 & \\ & 0 & \\ & & 0 \\ \end{bmatrix}.

Očito, treba nam prva derivacija funkcije f:

f'(t) = \cos t - e^{t}.

Sada možemo izračunati vrijednost f(A). Prema (4) i (5) imamo

\begin{align*} f(A)&= S^{-1} f\left(\begin{bmatrix} 0 & 1 & \\ & 0 & \\ & & 0 \\ \end{bmatrix}\right) S \\ &= \begin{bmatrix} -1 & -1/2 & -3 \\ 0 & 0 & 2 \\ 2 & 0 & 0 \end{bmatrix} \cdot \begin{bmatrix} f(0)& f'(0) \\ & f(0) \\ & & f(0) \\ \end{bmatrix} \cdot \begin{bmatrix} 0 & 0 & 1/2 \\ -2 & -3 & -1 \\ 0 & 1/2 & 0 \end{bmatrix} \\ &= S \cdot \begin{bmatrix} -1& 0 \\ & -1 \\ & & -1 \\ \end{bmatrix} \cdot S^{-1} = S^{-1} (-I) S = - S^{-1} S = -I. \end{align*}

2.2Hermiteov interpolacijski polinom

Kao i u prethodnom poglavlju, navest ćemo osnovne pojmove i rezultate, a zatim i definiciju matrične funkcije s pomoću Hermiteova interpolacijskog polinoma.

2.2.1Uvodni pojmovi i rezultati

Minimalni polinom matrice A jedinstveni je monični polinom \psi najmanjeg stupnja takav da je \psi(A) = 0. Nultočke minimalnog polinoma svojstvene su vrijednosti matrice i vrijedi

(6)

\psi(t) = \prod_{i=1}^{s} (t - \lambda_{i})^{n_{i}}.

Iz (6) je očito da je

\psi^{(j)}(\lambda_{i}) = 0, \quad j=0,\dots,n_{i}-1, \quad i=1,\dots,s,

odnosno, da su vrijednosti polinoma \psi i njegovih derivacija jednake nuli na cijelom spektru matrice. Ako je p polinom takav da je p(A) = 0 (naprimjer, karakteristični polinom), tada \psi dijeli p.

Nadalje, bilo koji polinom p definiran je na spektru od A, jer za polinome postoje sve derivacije u svim točkama. Posebno, to znači da postoje p^{(j)}(\lambda_{i}) za sve potrebne j i i.

Zanimljivo je primijetiti da su te vrijednosti dovoljne za jednoznačno određivanje p(A), tj.da vrijedi

Teorem 5. Neka su p i q polinomi i A\in \mathbb{C}^{n \times n}. Tada je p(A) = q(A) ako i samo ako je

p^{(j)}(\lambda_{i}) = q^{(j)}(\lambda_{i}), \quad j=0,\dots,n_{i}-1, \quad i=1,\dots,s.

Teorem je dokazan u Teoremu 1.3 u [1].

2.2.2Definicija matrične funkcije

Kao što smo vidjeli, lijepo svojstvo polinoma je da je matrica p(A) potpuno određena vrijednostima polinoma p (i njegovih prvih n_{i}-1 derivacija) na spektru od A. Prirodno generaliziramo to svojstvo na proizvoljnu funkciju f tako da definiramo f(A) na način da ta matrica bude potpuno određena samo vrijednostima funkcije f na spektru od A.

Definicija 6. (Matrična funkcija preko Hermiteova interpolacijskog polinoma) Neka je f definirana na spektru matrice A\in \mathbb{C}^{n \times n} i neka je

\psi(t) = \prod_{i=1}^{s} (t-\lambda_{i})^{n_{i}}

minimalni polinom od A. Tada je f(A) := p(A), gdje je p polinom stupnja manjeg ili jednakog \deg \psi = \sum_{i=1}^{s} n_{i} koji zadovoljava interpolacijske uvjete

(7)

p^{(j)}(\lambda_{i}) = f^{(j)}(\lambda_{i}), \quad j=0,\dots,n_{i}-1, \quad i=1,\dots,s.

Polinom p iz prethodne definicije je jedinstven i zove se Hermiteov interpolacijski polinom.

Ako je q bilo koji polinom (bilo kojeg stupnja) koji zadovoljava interpolacijske uvjete (7) iz definicije 6, tada se p iz definicije i q podudaraju na spektru od A, pa prema teoremu 5 slijedi da je q(A) = p(A) = f(A), tj. (7) jednoznačno određuje f(A).

Kao što vidimo, moguće je uzeti polinom q većeg stupnja nego što je nužno. Naprimjer, ako su poznate svojstvene vrijednosti, ali ne i Jordanova kanonska norma (dakle, pripadni indeksi svojstvenih vrijednosti), možemo uzeti n-1 derivacija funkcije za svaku svojstvenu vrijednost (umjesto njih n_{i}-1) i tako definirati potrebni polinom, te izračunati f(A) bez traženja indeksa svojstvenih vrijednosti.

Hermiteov interpolacijski polinom eksplicitno je dan Lagrange-Hermiteovom formulom:

(8)

p(t) = \sum_{i=1}^{s} \left[ \left( \sum_{j=0}^{n_{i}-1} \frac{\phi_{i}^{(j)}(\lambda_{i})}{j!}(t-\lambda_{i})^{j} \right) \prod_{\substack{j=1 \\ j \ne i}}^{n} (t - \lambda_{j})^{n_{j}} \right],

pri čemu je

\phi_{i}(t) = \displaystyle \frac{f(t)}{\displaystyle\prod_{\substack{j=1 \\ j \ne i}}^{n}(t-\lambda_{j})^{n_{j}}}.

Za matricu A kojoj su sve svojstvene vrijednosti međusobno različite (tj.s=n i n_{i}=1, za sve i), formula (8) prelazi u

p(t) = \sum_{i=1}^{n} f(\lambda_{i}) \prod_{\substack{j=1 \\ j \ne i}}^{n} \frac{t-\lambda_{j}}{\lambda_{i}-\lambda_{j}},

što je oblik Lagrangeova interpolacijskog polinoma.
Primijetimo da korištenjem definicije 6 lako možemo dobiti vrijednost f(J_{k}(\lambda_{k})) identičnu onoj iz (5).

2.3Ekvivalentnost dviju definicija

Prirodno se postavlja pitanje: koju od tih definicija odabrati? Srećom, kako je najavljeno u uvodu, dvije definicije su ekvivalentne, pa možemo birati onu koja nam u danom trenutku više odgovara.

Teorem 7. Funkcije f(A) definirane preko Jordanove kanonske forme i Hermiteova interpolacijskog polinoma međusobno su ekvivalentne.

Teorem je dokazan u Teoremu 1.12 u [1].

Svaka od dviju danih definicija za funkciju matrice f(A) ima svojih prednosti. Ako nam je poznata Jordanova kanonska forma matrice, iz definicije 2 relativno jednostavno dobivamo f(A).

S druge strane, definicija 6 eksplicitno kaže da je f(A) polinom u A. To znači da neovisno o izboru definicije funkcije f, funkciju matrice f(A) možemo izraziti kao polinom u A. Iz definicije 6 također slijedi:

(9)

Af(A) = f(A)A.

3Matrična funkcija predznaka

U ovom poglavlju primjenjujemo definicije iz prethodnog poglavlja na funkciju predznaka. Dat ćemo i kraći prikaz dvaju korištenih algoritama za računanje \mathop{\text{sign}} A te ćemo se osvrnuti na numeričko testiranje.

3.1Definicija \mathop{\text{sign}} A

Podsjetimo se: kao definiciju funkcije \mathop{\text{sign}} na skupu kompleksnih brojeva koristimo (3):

\mathop{\text{sign}}:\ \mathbb{C} \setminus \mathbb{I} \to \mathbb{R}, \quad \mathop{\text{sign}} z = \begin{cases} 1,& \Re z \gt 0, \\ -1,& \Re z \lt 0, \end{cases} \Bigg\rbrace = \mathop{\text{sign}}(\Re z).

Napomena 8.[Što s nulom?] Primijetimo da formula (3) ne definira vrijednosti funkcije predznaka na skupu imaginarnih brojeva. Razlog tomu je taj što nam za matričnu definiciju, u općenitom slučaju, trebaju derivacije funkcije. Ovako definirana funkcija je neprekidna (dapače, derivabilna), a svako dodefiniravanje na skupu imaginarnih brojeva bi, očito, uvelo prekide.

Iz (3) slijedi da matričnu funkciju predznaka možemo računati samo za matrice koje nemaju čisto imaginarne svojstvene vrijednosti, jer je funkcija predznaka definirana na njihovu spektru.

Lako se vidi da su sve derivacije funkcije \mathop{\text{sign}} jednake nuli, tj.

(10)

\text{sign}^{(k)}x = 0 \quad \text{za sve }x \in \mathbb{C} \setminus \mathbb{I}\text{}.

Iz (5) i (10) vidimo da se \mathop{\text{sign}} J_{k_{i}}(\lambda_{i}) računa podjednako jednostavno kao i za kompleksne brojeve:

(11)

\begin{aligned} \mathop{\text{sign}} J_{k_{i}}(\lambda_{i})&= \begin{bmatrix} \mathop{\text{sign}}(\lambda_{i}) \\ & \mathop{\text{sign}}(\lambda_{i}) \\ && \ddots \\ &&& \mathop{\text{sign}}(\lambda_{i}) \end{bmatrix} \\ &= \mathop{\text{sign}}(\lambda_{i}) \cdot I_{k_{i}} = \pm I_{k_{i}} \in \mathbb{C}^{k_{i} \times k_{i}}. \end{aligned}

Neka je A = S^{-1}JS Jordanova kanonska forma matrice A. Prema teoremu 1, Jordanove blokove možemo poredati po volji, pa ćemo ih grupirati tako da je J = \mathop{\text{diag}}(J_{1},J_{2}), pri čemu su Jordanovi blokovi koji odgovaraju svojstvenim vrijednostima u lijevoj poluravnini (realni dio im je negativan) grupirani u J_{1}, a oni koji odgovaraju svojstvenim vrijednostima u desnoj poluravnini (realni dio im je pozitivan) grupirani u J_{2}.

Iz definicije 2 i formule (11) slijedi da je

(12)

\mathop{\text{sign}} A = S^{-1} \mathop{\text{sign}}(J) S = S^{-1} \mathop{\text{sign}}(\mathop{\text{diag}}(J_{1}, J_{2}) S = S^{-1} \mathop{\text{diag}}(-I_{p}, I_{q}) S,

pri čemu je p ukupna algebarska kratnost svih svojstvenih vrijednosti u lijevoj poluravnini (red matrice J_{1}), a q ukupna algebarska kratnost svih svojstvenih vrijednosti u desnoj poluravnini (red matrice J_{2}).

Iz (12) vidimo da je matrična funkcija \mathop{\text{sign}} involutorna:

(13)

\begin{aligned} (\mathop{\text{sign}} A)^{2}&= (S^{-1} \mathop{\text{diag}}(-I_{p}, I_{q}) S)^{2} = S^{-1} \mathop{\text{diag}}(-I_{p}, I_{q}) S \cdot S^{-1} \mathop{\text{diag}}(-I_{p}, I_{q}) S \\ &= S^{-1} (\mathop{\text{diag}}(-I_{p}, I_{q}))^{2} S = S^{-1} (\mathop{\text{diag}}((-I_{p})^{2}, I_{q}^{2})) S \\ &= S^{-1}\,I\,S = I. \end{aligned}

Napomena 9. Ako su sve svojstvene vrijednosti matrice A u lijevoj poluravnini (tj.p = n i q = 0), onda vrijedi: \mathop{\text{sign}} A = -I. Analogno, ako su sve svojstvene vrijednosti matrice A u desnoj poluravnini (tj.p = 0 i q = n), onda vrijedi: \mathop{\text{sign}} A = I.

3.2Računanje \mathop{\text{sign}} A

Teoretska definicija matričnih funkcija (i, posebno, matrične funkcije predznaka) omogućuje dobro razumijevanje pojedinih funkcija, no ostaje problem konkretnih algoritama za izračunavanje njihovih vrijednosti. U ovom poglavlju razmotrit ćemo Schurovu metodu računanja vrijednosti \mathop{\text{sign}} A za zadanu matricu A, te dati kraći opis Newtonove metode za računanje iste vrijednosti.

3.2.1Schurova metoda

Metoda koju ćemo prikazati bazirana je na Schurovoj dekompoziciji koja kaže da je svaka kvadratna matrica unitarno slična nekoj gornjotrokutastoj matrici. Za matricu U\in \mathbb{C}^{n \times n} kažemo da je unitarna ako je U^{*}U = UU^{*} = I. Očito, takva matrica je regularna i vrijedi: U^{-1} = U^{*}.

Teorem 10. (Schurova dekompozicija) Neka je A\in \mathbb{C}^{n \times n} kompleksna matrica s (ne nužno različitim) svojstvenim vrijednostima \lambda_{1}, \lambda_{2},\dots, \lambda_{n}. Tada postoje unitarna matrica U\in \mathbb{C}^{n \times n} i gornjotrokutasta matrica T \in \mathbb{C}^{n \times n} takve da je

(14)

A = U^{*}TU \quad \text{i} \quad t_{ii} = \lambda_{i}.

Ako je A = S^{-1}BS i ako je B = Z^{-1}JZ Jordanova kanonska forma matrice B, onda za vrijednosti funkcije f(A) i f(B), korištenjem relacije (4), slijedi:

\begin{align*} f(A)&= f(S^{-1}BS) = f(S^{-1}Z^{-1}JZS) = f((ZS)^{-1}J(ZS)) \\ &= (ZS)^{-1}f(J)(ZS) = S^{-1}Z^{-1}f(J)ZS = S^{-1}(Z^{-1}f(J)Z)S \\ &= S^{-1}f(Z^{-1}JZ)S = S^{-1}f(B)S. \end{align*}

Zbog U^{*} = U^{-1}, prva relacija u (14) je sličnost matrica, pa prema upravo pokazanoj povezanosti funkcijske vrijednosti sličnih matrica, za matrice A, U i T iz teorema 10 vrijedi:

(15)

f(A) = U^{*}f(T)U.

Drugim riječima, ako je poznata Schurova dekompozicija matrice, računanje f(A) svodi se na računanje vrijednosti funkcije gornjotrokutaste matrice. Može se pokazati da vrijedi sljedeći teorem.

Teorem 11. Ako je T\in \mathbb{C}^{n \times n} gornjotrokutasta matrica i f funkcija definirana na spektru od T, tada je F = f(T) također gornjotrokutasta matrica s dijagonalnim elementima f_{ii} = f(t_{ii}).

Teorem je dokazan u Teoremu 1.13 f) [1].

Primjenom (15) na računanje \mathop{\text{sign}} A, vidimo da se ono svodi na računanje \mathop{\text{sign}} T, jer je

(16)

\mathop{\text{sign}} A = U^{*}(\mathop{\text{sign}} T)U.

Prema teoremu 11, vidimo da je matrica S := \mathop{\text{sign}} T gornjotrokutasta i za njene dijagonalne elemente vrijedi:

s_{ii} = \mathop{\text{sign}}(t_{ii}) = \pm1, \ \text{ za svaki }i\text{}.

Ostaje izračunati vrijednosti onih elemenata matrice S koji se nalaze iznad dijagonale — tražimo vrijednosti s_{ij} za i\lt j.

Te vrijednosti izračunat ćemo korištenjem involutornosti \mathop{\text{sign}} A (formula (13)) i TS = ST, što slijedi direktno iz (9). U nastavku ćemo se koristiti Kroneckerovom deltom:

\delta_{ij} = \begin{cases} 1,& i = j, \\ 0,& i \ne j. \end{cases}

Pogledajmo prvo elemente matrice S^{2}, za koju — zbog involutornosti funkcije \mathop{\text{sign}} — vrijedi S^{2} = I:

\delta_{ij} = [S^{2}]_{ij} = \begin{cases} 0,& i \gt j, \\ s_{jj}s_{jj} = s_{jj}^{2},& i = j, \\ \displaystyle \sum_{k=i}^{j} s_{ik}s_{kj},& i \lt j. \end{cases}

Za i \geqslant j nismo saznali ništa novo, ali za i \lt j imamo:

\sum_{k=i}^{j} s_{ik}s_{kj} = s_{ii}s_{ij} + s_{i,i+1}s_{i+1,j} + \dots + s_{ij}s_{jj} = \delta_{ij} = 0.

Element s_{ij} koji želimo izračunati pojavljuje se u prvom i zadnjem članu sume. Ako je s_{ii} + s_{jj} \ne 0, lako vidimo da je

(17)

s_{ij} = \frac{\displaystyle -\sum_{k=i+1}^{j-1}s_{ik}s_{kj}}{s_{ii}+s_{jj}}.

Primijetimo da nam za računanje s_{ij} po formuli (17) — ako je s_{ii} + s_{jj} \ne 0 — trebaju elementi matrice S koji se nalaze između glavne dijagonale i samog elementa s_{ij} (tj.u istom retku lijevo i istom stupcu ispod njega). To znači da elemente matrice S možemo računati redom, bilo paralelno s glavnom dijagonalom, bilo po stupcima. U nastavku ćemo prikazati računanje po stupcima.

Preciznije, prvo računamo elemente s_{ii} na glavnoj dijagonali, zatim preostali element drugog stupca s_{12}, zatim preostala dva elementa trećeg stupca s_{23} i s_{13} itd. Važno je da elemente pojedinog stupca računamo od najdonjeg prema najgornjem, dakle

s_{j-1,j},\ s_{j-2,j}, \ \dots, \ s_{1,j}

jer nam za računanje pojedinog elementa trebaju biti izračunati svi elementi ispod njega. Iz formule (17) vidimo i da je za elemente s_{i,i+1} suma u toj formuli prazna (dakle, jednaka nuli), pa vrijedi:

s_{ii} + s_{i+1,i+1} \ne 0 \quad \Rightarrow \quad s_{ij} = 0.

Ako se u nekom koraku dogodi s_{ii} + s_{jj} = 0, onda se ne možemo koristiti formulom (17). Primijetimo još da će se to sigurno dogoditi u svakoj matrici koja nema sve svojstvene vrijednosti u istoj (lijevoj ili desnoj) poluravnini, jer je s_{ii} = \mathop{\text{sign}}(\lambda_{i}).

U tom slučaju koristimo se jednakošću TS = ST. Budući da su i T i S gornjotrokutaste matrice (T prema teoremu 10, a S prema teoremu 11), onda su i njihovi produkti TS i ST također gornjotrokutaste matrice. Dijagonalni elementi tih produkata su t_{ii}s_{ii}, dok su im poddijagonalni elementi jednaki nuli.

Zanimljivi za promatranje su naddijagonalni elementi s_{ij}, za i \lt j. Prema definiciji produkta matrica, iz TS = ST slijedi da je

\sum_{k=1}^{j} t_{ik}s_{kj} = \sum_{k=1}^{j} s_{ik}t_{kj}.

Jednostavnim raspisom dobijemo:

(18)

t_{ii}-t_{jj} \ne 0 \quad \Rightarrow \quad s_{ij} = \frac{t_{ij}(s_{ii}-s_{jj}) + \displaystyle \sum_{k=i+1}^{j-1}(s_{ik}t_{kj}-t_{ik}s_{kj})}{t_{ii}-t_{jj}}.

Očito, sve napomene vezane uz formulu (17) vrijede i za formulu (18). Redoslijed računanja jednak je kao i za formulu (17), a suma u formuli prazna je za sve elemente neposredno iznad dijagonale.
Primijetimo da vrijedi:

\begin{align*} t_{ii}-t_{jj} = 0&\Longrightarrow t_{ii} = t_{jj} \Longrightarrow s_{ii} = \mathop{\text{sign}} t_{ii} = \mathop{\text{sign}} t_{jj} = s_{jj} \in \lbrace -1,1\rbrace \\ &\Longrightarrow s_{ii} + s_{jj} \in \lbrace -2,2\rbrace \\ &\Longrightarrow s_{ii} + s_{jj} \ne 0. \end{align*}

Drugim riječima, ako formula (18) nije primjenjiva, onda se sigurno možemo koristiti formulom (17). To znači da su pokriveni svi slučajevi, pa su formule (17) i (18) dovoljne za računanje vrijednosti funkcije predznaka za gornjotrokutaste matrice. Naravno, ako možemo birati koju formulu upotrijebiti, očito ćemo – radi jednostavnosti – upotrijebiti (17).
Osim funkcije predznaka, s pomoću formule (18) mogu se računati vrijednosti bilo koje funkcije na gornjotrokutastim matricama, uz uvjet da je t_{ii} \ne t_{jj}, za sve različite i i j (tj.da matrica nema višestrukih svojstvenih vrijednosti), jer se a formula nigdje ne koristi ničim specifičnim za funkciju predznaka. Za razliku od nje, formula (17) se koristi idempotentnošću funkcije predznaka, što je svojstvo koje nemaju sve funkcije.

Pogledajmo sada kako izgleda algoritam kojim se računa \mathop{\text{sign}} A proizvoljne kvadratne matrice A\in \mathbb{C}^{n \times n}.

Algoritam 12.[Schurova metoda] Za A\in \mathbb{C}^{n \times n} koja nema čisto imaginarnih svojstvenih vrijednosti, algoritam računa F = \mathop{\text{sign}} A korištenjem Schurove dekompozicije stupac po stupac.

(1)	Odrediti (kompleksnu) Schurovu dekompoziciju A = U^{*}TU.
(2)	Za sve i = 1,2,\dots,n, izračunati s_{ii} = \mathop{\text{sign}} t_{ii}.
(3)	Za sve stupce j=2,3,\dots,n gledati retke i=j-1,j-2,\dots,1 i računati s_{ij} = \begin{cases} \displaystyle \frac{\displaystyle -\sum_{k=i+1}^{j-1}s_{ik}s_{kj}}{s_{ii}+s_{jj}},& s_{ii}+s_{jj} \ne 0, \\ \displaystyle \frac{t_{ij}(s_{ii}-s_{jj}) + \displaystyle \sum_{k=i+1}^{j-1}(s_{ik}t_{kj}-t_{ik}s_{kj})}{t_{ii}-t_{jj}},& \text{inače}. \end{cases}
(4)	Izračunati F = U^{*}SU.

Složenost prikazanog algoritma je kubna, a broj potrebnih računskih operacija nad kompleksnim brojevima iznosi približno \frac{86}{3}n^{3}. Detaljni raspis ove tvrdnje dostupan je na stranicama 22.–25. u [6].

3.2.2Newtonova metoda

Kod Schurove metode iterativnost je donekle sakrivena, u smislu da sama metoda nije iterativna, ali Schurova dekompozicija jest. Za razliku od nje, Newtonova metoda je očigledno iterativna.

Osnovna ideja je računanje prvih nekoliko elemenata niza matrica A_{n}, pri čemu taj niz definiramo na način da brzo konvergira traženoj vrijednosti f(A). Uz računanje elemenata niza, računa se i procijenjeno odstupanje od tražene vrijednosti koje se koristi kao kriterij zaustavljanja računanja.

Metoda se zove Newtonova, jer je, u osnovi, riječ o traženju nultočke funkcije Newtonovom metodom. Naprimjer, za funkciju predznaka znamo da je involutorna, tj.da za F = \mathop{\text{sign}} A vrijedi: F^{2} = I, odnosno F^{2} - I = 0. Da je riječ o brojevima, tražili bismo nultočku funkcije f(x) := x^{2} - 1.

Newtonova metoda za skalarne jednadžbe f(x) = 0 kaže da niz elemenata x_{k} za koje vrijedi

x_{k+1} = x_{k} - \frac{x_{k}^{2}-1}{2x_{k}} = \frac{1}{2}\left(x_{k}+\frac{1}{x_{k}}\right) = \frac{1}{2}(x_{k}+x_{k}^{-1}),

uz neke dodatne uvjete, konvergira traženom rješenju x_{0} jednadžbe f(x)=0. Prethodnu formulu prirodno prepisujemo u matrični oblik:

A_{k+1} = \frac{1}{2}(A_{k}+A_{k}^{-1}).

Algoritam očito treba računati prvih nekoliko članova niza (A_{k}), dok ne dođe do onoga koji je zadovoljavajuće blizu rješenju. Postavlja se pitanje koji je to član, tj.koji kriterij zaustavljanja računanja iteracija primijeniti?

U praksi se koriste različiti kriteriji zaustavljanja, bazirani na matričnim normama. Uobičajeno je promatrati relativnu razliku dviju susjednih iteracija,

\delta_{k+1} := \frac{\Vert A_{k+1}-A_{k}\Vert }{\Vert A_{k+1}\Vert },

pri čemu je \Vert \cdot\Vert neka matrična norma.

U našem slučaju, želimo F := \lim_{k \to \infty} A_{k} aproksimirati nekim A_{k+1} tako da relativna greška za A_{k+1} bude manja od unaprijed zadanog \varepsilon \gt 0. No, limes F nam nije poznat. Zato možemo promatrati \delta_{k+1} koji daje informaciju „koliko se rezultat mijenja između dviju iteracija”. Kad promjena postane dovoljno mala, smatramo da je rješenje zadovoljavajuće.

Budući da kriterij zaustavljanja treba računati u svakom koraku, važno je da to računanje bude što brže. Pogodne matrične norme za takav račun su, naprimjer,

(1)	1-norma: \Vert A\Vert _{1} := \displaystyle \max_{1 \leq j \leq n} \sum_{i=1}^{n} |a_{ij}|,
(2)	Frobeniusova ili F-norma: \Vert A\Vert _{\text{}}{F} := \displaystyle \sqrt{\sum_{i=1}^{n} \sum_{j=1}^{n} |a_{ij}|^{2}},
(3)	\infty-norma: \Vert A\Vert _{\infty} := \displaystyle \max_{1 \leq i \leq n} \sum_{j=1}^{n} |a_{ij}|.

Umjesto gledanja općenite relativne greške, za kriterij zaustavljanja iskoristit ćemo svojstvo kvadratne konvergencije Newtonovih iteracija

\Vert A_{k+1} - F\Vert \leqslant c \Vert A_{k} - F\Vert ^{2},

za proizvoljnu matričnu normu (uzet ćemo Frobeniusovu), jer ona daje bolju ocjenu greške i, samim tim, bolji kriterij zaustavljanja.
Jednostavnim raspisom dobivamo kriterij zaustavljanja baziran na Frobeniusovoj normi:

(19)

\mathop{\text{s}top}(F, k, \varepsilon) \equiv \left( \Vert F_{k+1}-F_{k}\Vert _{\text{}}{F} \leqslant \sqrt{\frac{\varepsilon \Vert F_{k+1}\Vert _{\text{}}{F}}{\Vert F_{k}^{-1}\Vert _{\text{}}{F}}} \right).

Kod iterativnih algoritama dobro je osigurati se od divergencije. Iako algoritam u teoriji uvijek konvergira, zbog grešaka konačne računalne aritmetike može se dogoditi da program u nekim slučajevima divergira. Kako program niti u tom slučaju ne bi “zapeo” u beskonačnoj petlji, nužno je unaprijed zadati gornju ogradu za broj dopuštenih iteracija. Primjere takve divergencije navest ćemo u poglavlju 3.3.

Sada možemo prikazati i algoritam za računanje vrijednosti \mathop{\text{sign}} A Newtonovom metodom, za zadanu matricu A.

Algoritam 13.[Newtonova metoda] Za A\in \mathbb{C}^{n \times n} koja nema čisto imaginarnih svojstvenih vrijednosti i za zadanu točnost \varepsilon \gt 0, te najveći dopušteni broj iteracija, algoritam računa F = \mathop{\text{sign}} A.

(1)	F_{0} = A i k = 0.
(2)	Dok je k \lt \text{max_ {i}teracija}[1.] (c) G_{k} = F_{k}^{-1} (d) F_{k+1} = (F_{k}+G_{k})/2 (e) Ako je \mathop{\text{s}top}(F, k, \varepsilon) = \verbISTINA, prekini petlju; inače povećaj k.
(6)	F = F_{k+1}.

Postavlja se logično pitanje: zašto se koristiti očigledno kompliciranijom Newtonovom metodom umjesto Schurove? Razlog leži u brzini metode, kao i činjenici da možemo kontrolirati postignutu točnost rješenja.

Newtonova metoda će u k koraka obaviti otprilike 2kn^{3} računskih operacija, što znači da će tek s 15 ili više koraka iteracije postati sporija od Schurove metode. U praksi, iterativni algoritmi obično brzo konvergiraju i 15 koraka uopće nije malo, a broj potrebnih iteracija može se dodatno smanjiti uvođenjem skaliranja u algoritam (v.poglavlje 5 u 3.3).

3.3Numeričko ispitivanje

Numeričko ispitivanje, prikazano u poglavlju 5 vsego-diplNat}, napravljeno je na matricama reda n \leqslant 1000 bez čisto imaginarnih svojstvenih vrijednosti, korištenjem biblioteke MKL3 za Linux, verzija 8.0.2, na računalu

DELL PowerEdge 2800
Procesori: 2 x Intel Xeon, 3.0GHz/1MB L2 Cache
RAM: 2GB 400MHz ECC DDR2
HDD: 4 x 146GB Hot-Plug Ultra 320 SCSI RAID 5, 256MB cache
NIC: Dual Intel Gigabit
OS: Debian Lenny (5.0.1).

Testiranje je podijeljeno u dvije grupe. Prvu čine tzv. normalne matrice4. Za potrebe testiranja, bitno je samo da su to matrice čija je Schurova dekompozicija oblika U^{*}DU, pri čemu je D dijagonalna matrica. Tada i Schurova i Newtonova metoda daju rezultate visoke točnosti, s time da je Newtonova donekle sporija (za matrice reda oko 700 traje otprilike dvostruko dulje od Schurove metode i ta razlika polagano raste s redom matrice).

Druga grupa promatranih matrica su one čija je Schurova dekompozicija oblika U^{*}TU, pri čemu je T gornjotrokutasta nedijagonalna matrica (tj.postoje i i j takvi da je i\lt j i t_{ij} \ne 0). Za takve matrice Schurova metoda opet daje zadovoljavajuće rezultate, dok Newtonova metoda gotovo uvijek divergira.

Veliki problem Newtonove metode je što u svakom koraku računa matrični inverz. Za neke (tzv. loše uvjetovane) matrice inverz izračunat na računalu može biti jako daleko od stvarnog inverza. Ispitivanje u Mathematici5 pokazalo je da do divergencije u računanju \mathop{\text{sign}} A dolazi upravo zbog velikih grešaka u računanju inverza loše uvjetovanih matrica.

Slika 1: Divergencija Newtonove metode za ne-normalne matrice

Na slici 1 prikazano je računanje \mathop{\text{sign}} A za jednu takvu slučajnu matricu reda 50 u prvih 1000 koraka. Na ordinati je prikazana “udaljenost” (izračunata s pomoću Frobeniusove norme) od referentnog rješenja (izračunato direktno po definiciji, tj.preko Jordanove kanonske forme). Točni brojevi nisu toliko bitni koliko očita činjenica da algoritam ne konvergira (barem ne u prvih 1000 koraka). Na slici nije prikazano, ali pokazalo se da metoda ne konvergira niti u 10^{4} koraka. Za usporedbu: kad algoritam konvergira, onda mu je prosječan broj koraka oko 15 za matricu reda 50.

Zanimljivo je ove metode usporediti i na jednom “namještenom” primjeru. Riječ je o tzv. Hilbertovim matricama, tj. matricama oblika

H_{n} := [ h_{ij} ], \quad h_{ij} = \frac{1}{i+j-1}, \quad i,j \in \lbrace 1,2,\dots,n\rbrace .

Očito, Hilbertove su matrice hermitske, što znači i da su normalne. Spektar prikazane matrice H_{10} je

\begin{align*} \sigma(H_{10}) = \lbrace &1.75192, 0.34293, 0.0357418, 0.00253089, 0.00012875, \\ &4.72969\cdot10^{-6}, 1.22897\cdot10^{-7}, 2.14744\cdot10^{-9}, \\ &2.26674\cdot10^{-11}, 1.09312\cdot10^{-13}\rbrace . \end{align*}

Kao što vidimo, sve svojstvene vrijednosti matrice H_{10} su pozitivne (iako jako blizu nuli), što, prema napomeni 9, znači da je \mathop{\text{sign}} H_{10} = I.

Schurova metoda rezultat daje iznimno brzo i s velikom točnošću. Konkretno, za H_{10} rezultat je dobiven za samo 0.004 sekunde, uz točnost \Vert I - \mathop{\text{sign}} H_{10}\Vert _{\text{}}{F} = 1.746977\cdot10^{-15}.

Detaljni rezultati testiranja Newtonove metode na matrici H_{10} dani su u tablici 3 u [6]. Metoda je također brza, ali je broj iteracija izrazito velik. Za točnost reda veličine 10^{-10}, potrebno je čak 47 iteracija. Usporedbe radi, za slučajne normalne matrice reda 10, prosječni broj iteracija je približno 11.4.

Zbog svega opisanog, za računanje funkcije predznaka Schurova je metoda prikladnija od osnovne Newtonove. S obzirom na probleme oko računanja matričnog inverza, za očekivati je da niti skaliranje ne bi bitno popravilo rezultate dobivene Newtonovom metodom.

4Zaključak

U ovom radu dali smo motivaciju i uvod u teoriju matričnih funkcija. Prikazane su dvije ekvivalentne definicije, svaka sa svojim prednostima, te ispitana osnovna svojstva matričnih funkcija. U nastavku smo detaljnije promotrili matričnu funkciju predznaka, a zatim je detaljno opisana Schurova metoda za računanje vrijednosti \mathop{\text{sign}} A. Opisali smo i Newtonovu metodu, te dali kraću usporedbu dviju metoda.

Vidjeli smo da prelazak s “običnih” funkcija (definiranih za realne ili kompleksne brojeve) na matrične funkcije nije trivijalan. Dodatno, kod same funkcije predznaka morali smo se odreći dijela domene (matricâ s čisto imaginarnim svojstvenim vrijednostima) zbog dodatnog zahtjeva definicije matrične funkcije na derivabilnost početne funkcije u točkama u kojima matrica ima svojstvene vrijednosti.

Kao i sva prirodna proširenja domena funkcija, matrična proširenja funkcija koncentriraju se na dobra svojstva funkcija koja zatim pokušavaju zadržati i u svojim matričnim varijantama. Kao primjer toga vidjeli smo da kod funkcije predznaka čuvamo idempotentnost, zahvaljujući kojoj smo dobili i bitno jednostavniju formulu (17) za računanje vrijednosti funkcije.

Matrična funkcija predznaka posebno je jednostavna jer su sve derivacije funkcije \mathop{\text{sign}} jednake nuli. No, unatoč toj jednostavnosti, ona ima mnoge primjene. Naprimjer, s pomoću vrijednosti \mathop{\text{sign}} A, prema (12), vrlo jednostavno možemo odrediti broj svojstvenih vrijednosti matrice A u lijevoj i desnoj poluravnini. Jednostavnom translacijom A_{c} := A - cI, s pomoću vrijednosti \mathop{\text{sign}} A_{c} možemo lako odrediti broj svojstvenih vrijednosti matrice A čiji je realni dio manji, odnosno veći od c. Pogodnim odabirom parametra c možemo dobiti broj svojstvenih vrijednosti matrice A čiji je realni dio u nekom intervalu, što je u praksi korisno za velike matrice, kada nas zanima samo dio spektra.

Rad se bazira na diplomskom radu N.Strabić [6] i u njemu se mogu naći detaljno raspisane sve korištene tvrdnje. Zainteresirani čitatelji u [1] mogu naći puno detaljnije razrađenu teoriju matričnih funkcija, uz detaljan pregled drugih matričnih funkcija (poput matričnih potencija, korijena i logaritama), ocjene greške, racionalne aproksimacije, matrične dekompozicije vezane uz računanje matričnih funkcija i sl.

5Zahvale

Ovim putem zahvaljujemo prof.dr.sc.Saši Singeru na velikoj pomoći s debuggiranjem i testiranjem programa i doc.dr.sc.Lavoslavu Čakloviću na korištenju računala na kojem su testiranja provedena.

Također, Ines Šimičić i prof.dr.sc.Sanji Singer zahvaljujemo na generalnim smjernicama oko same radnje na temelju koje je nastao ovaj rad.

Bibliografija

[1]	N.J. Higham. Functions of matrices. Theory and computation. Philadelphia, PA: Society for Industrial and Applied Mathematics (SIAM). 425p., 2008
[2]	R.A. Horn and C.R. Johnson. Matrix Analysis. Cambridge University Press, Cambridge, 1985
[3]	K.Horvatić. Linearna algebra, II. dio. Matematički odjel PMF–a Sveučilišta u Zagrebu, Zagreb, 1995
[4]	K.Horvatić. Linearna algebra, III. dio. Matematički odjel PMF–a Sveučilišta u Zagrebu, Zagreb, 1995
[5]	S.Kurepa. Konačno dimenzionalni vektorski prostori i primjene. Tehnička knjiga, Zagreb, 1967
[6]	N.Strabić. Matrična funkcija predznaka. PMF – Matematički odjel, Sveučilište u Zagrebu, lipanj 2009 diplomski rad, http://viveka.math.hr/~vsego/nstrabic-dipl.

---

1Za matricu T = [ t_{ij} ]\in \mathbb{C}^{n \times n} kažemo da je gornjetrokutasta ako za sve i\gt j vrijedi da je t_{ij}=0, tj.ako su svi elementi ispod glavne dijagonale jednaki nuli, dok na dijagonali i iznad nje mogu biti proizvoljni elementi (uključujući i nule).

2Naprimjer, ako je S matrica sličnosti, tada je to i svaka matrica \alpha S, za \alpha \neq 0.

3Intel® Math Kernel Library

4Za kvadratnu matricu N kažemo da je normalna ako vrijedi N^{*}N=NN^{*}.

5©Wolfram Research, Inc.Mathematica 5.1.0.0 za Linux