Nieregularne LZNK

Jeśli rząd r macierzy A jest mniejszy zarówno od liczby kolumn n, jak i od liczby wierszy m, to liniowe zadanie najmniejszych kwadratów dla układu Ax = b jest nieregularne. Zbiór rozwiązań takiego zadania jest nieskończony; jest on warstwą n − r-wymiarową (przestrzeni Rⁿ), której elementami są takie wektory x, że wektor y^⋆= Axjest rzutem prostopadłym wektora b na podprzestrzeń rozpiętą

przez kolumny macierzy A (tj. residuum, b − y^⋆, jest wektorem prostopadłym do tej podprzestrzeni). Dokładnie jeden element tej warstwy ma najmniejszą normę drugą; co więcej, dla dowolnego wektora ^x ∈ Rⁿistnieje dokładnie jeden

element x^⋆tej warstwy, taki że norma druga różnicy x^⋆− ^xjest najmniejsza.

Rozwiązanie NLZNK zwykle polega na znalezieniu tego wektora x^⋆. Rozumowanie podobne do przeprowadzonego wcześniej dla DLZNK uzasadnia stwierdzenie, że wektor x^⋆− ^xjest kombinacją liniową transponowanych wierszy macierzy A.

NLZNK są trudne do numerycznego rozwiązania. Jest tak dlatego, że rozwiązanie zadania zależy od danych w sposób paskudnie nieciągły. NLZNK jest szczególnie trudne, jeśli nie znamy rzędu macierzy A i dopiero mamy go na podstawie obliczeń numerycznych ustalić.

Najodporniejsze numeryczne algorytmy rozwiązywania NLZNK korzystają z rozkładu względem wartości szczególnych (ang. singular value decomposition, SVD) macierzy A. Dowodzi się, że dla dowolnej macierzy A ∈ R^m,nistnieją macierze ortogonalne U ∈ R^m,mi V ∈ R^n,noraz macierz diagonalna Σ ∈ R^m,n, takie że A = UΣV^T. Współczynniki diagonalne macierzy Σ, σ1, . . . , σl, gdzie l =min{m, n}, nazywają się wartościami szczególnymi macierzy A i są nieujemne.

Rozkład w ogólności nie jest jednoznaczny, ale same wartości szczególne i liczby ich wystąpień (krotności) są określone przez macierz A jednoznacznie. Zwykle rozkładu dokonuje się w taki sposób, aby wartości szczególne były uporządkowane nierosnąco na diagonali macierzy Σ: σ1>σ2>· · · > σ^r> σr+1=· · · = σ^l= 0.

Liczba niezerowych wartości szczególnych jest rzędem macierzy A.

Wyznaczanie rozkładu SVD wiąże się z rozwiązywaniem algebraicznego

zagadnienia własnego i dla macierzy o więcej niż czterech wierszach i kolumnach może być dokonane tylko jakąś metodą iteracyjną. Opis algorytmu Goluba, który dokonuje rozkładu (otrzymując reprezentacje macierzy U i V w postaci ciągu wektorów odbić Householdera i tzw. obrotów Givensa), pominiemy (kto będzie potrzebował, ten go znajdzie). Natomiast przyjrzyjmy się zastosowaniu tego rozkładu do rozwiązania zadania.

NLZNK dla układu równań Ax = b i wektora ^x można zastąpić zadaniem równoważnym dla układu równań Σy = d, gdzie d = U^Tb, i wektora ^y = V^T^x (po rozwiązaniu tego zadania możemy obliczyć x = Vy). Załóżmy dla uproszczenia, że ^x = 0, czyli ^y = 0. Wtedy rozwiązaniem NLZNK dla układu

Σy = djest wektor o współrzędnych

yi=

 d_i/σ_i dla i 6 r, 0 dla i > r.

Mamy stąd wyjaśnienie trudności zadania: niewielkie zaburzenie macierzy A może spowodować pewne niegroźne zmiany macierzy U i V, oraz zaburzenie macierzy Σ:

jeśli dowolna zerowa wartość szczególna zmieni się na niezerową (czyli skutkiem zaburzenia będzie zwiększenie rzędu macierzy A) i di6= 0, to trzeba będzie przyjąć y_i= d_i/σ_i, zamiast zera, dla pewnego i > r. Tak więc, im mniej zaburzymy macierz A (w sposób zmieniający σi), tym większa będzie zmiana wyniku.

Jeśli znamy rząd r macierzy A, to po znalezieniu rozkładu SVD możemy zamienić na zera obliczone numerycznie wartości szczególne σidla i > r — obliczone wartości niezerowe są skutkiem błędów zaokrągleń i aproksymacji popełnionych podczas rozkładania. Jeśli rzędu nie znamy, to możemy przyjąć pewien próg (zależny od oszacowania błędów) i zamienić na zera znalezione wartości szczególne mniejsze od tego progu; to postępowanie nazywa się regularyzacją dyskretną.

Inne podejście to tzw. regularyzacja ciągła — do wszystkich wartości szczególnych dodajemy pewną liczbę s > 0, otrzymując zadanie z macierzą pełnego rzędu, tj. RLZNK, jeśli m > n, układ równań z macierzą kwadratową nieosobliwą, jeśli n = m, albo DLZNK, jeśli m < n. Wybór metody regularyzacji zależy od zastosowania.

Zadania i problemy

1. Niech

A =









 1 1 1 ε 0 0 0 ε 0 0 0 ε









 , b =









 3 ε ε ε









 .

Jaki będzie skutek użycia arytmetyki pojedynczej precyzji do rozwiązania LZNK dla układu Ax = b za pomocą algorytmu równań normalnych, jeśli |ε| < 10⁻⁴? 2. Tabela zawiera wyniki pomiarów (z błędami) wartości pewnej funkcji:

x −2 −1 1 2

f(x) −1 5 2 14

Przy założeniu, że funkcja ma postać f(x) = a0+ a₁(x²− 4), znajdź liczby a0i a1

najlepiej pasujące do wyników tych pomiarów. Postaw i rozwiąż w tym celu LZNK, używając algorytmu równań normalnych oraz odbić Householdera.

3. Pseudoodwrotność macierzy A ∈ R^m,njest to macierz A⁺∈ R^n,m, taka że macierz AA⁺jest macierzą rzutu ortogonalnego przestrzeni R^mna podprzestrzeń rozpiętą przez kolumny macierzy A, zaś macierz A⁺Ajest macierzą rzutu ortogonalnego na podprzestrzeń rozpiętą przez transponowane wiersze macierzy A.

Udowodnij, że

a) jeśli macierz A jest kwadratowa nieosobliwa, to A⁺= A⁻¹, b) jeśli macierz A jest kolumnowo regularna, to A⁺= (A^TA)⁻¹A^T, c) jeśli macierz A jest wierszowo regularna, to A⁺= A^T(AA^T)⁻¹.

Zbadaj związek pseudoodwrotności z liniowymi zadaniami najmniejszych kwadratów.

Wskazówka: znajdź wzory opisujące rozwiązania odpowiednich układów równań normalnych.

Uwaga: Pseudoodwrotność macierzy A czasem definiuje się jako macierz A⁺ spełniającą tzw. warunki Penrose’a:

A⁺A = (A⁺A)^T, AA⁺= (AA⁺)^T, AA⁺A = A, A⁺AA⁺= A⁺. Dowodzi się, że dla dowolnej macierzy A istnieje dokładnie jedna macierz A⁺ spełniająca te warunki. Ta definicja jest równoważna definicji podanej wcześniej w tym zadaniu. Ponadto, znając rozkład SVD macierzy A, tj. macierze

ortogonalne U, V oraz macierz diagonalną Σ, takie że A = UΣV^T, możemy napisać A⁺= VΣ⁺U^T, gdzie Σ⁺jest macierzą diagonalną ze współczynnikami

1/σ1, . . . , 1/σr, 0, . . . , 0na diagonali.

4. LZNK z więzami. Dany układ równań liniowych jest podzielony na dwa podukłady, z macierzami B ∈ R^m,ni C ∈ R^k,n(takimi, że k < n < m + k):

 Bx = d, Cx = e,

przy czym macierz A = ^B_C

jest kolumnowo regularna, a macierz C jest wierszowo regularna. Zadanie polega na znalezieniu wektora x^⋆, takiego że Cx^⋆= eoraz norma druga wektora residuum pierwszego podukładu jest najmniejsza.

Algorytm zamiany zmiennych: Znajdujemy rozkład macierzy C na czynniki L i Q^T, odpowiednio trójkątny dolny i ortogonalny. Dokonujemy zamiany

zmiennych, wprowadzając nowy wektor niewiadomy, y = Q^Tx. Macierz L dzielimy na bloki, L = [L1, L₂], z których pierwszy jest nieosobliwą macierzą trójkątną dolną, a drugi jest macierzą zerową. Macierz Q dzielimy na bloki Q1i Q2. W ten sposób drugi podukład zastępujemy układem równoważnym L1y₁= e, z którego wyznaczymy y1. Po zamianie zmiennych w pierwszym podukładzie otrzymujemy układ równań liniowych

BQy = d, czyli BQ1y1+ BQ₂y2= d, czyli BQ2y2= d − BQ₁y1. w ogólności sprzeczny. Dla tego układu rozwiązujemy regularne liniowe zadanie najmniejszych kwadratów, po czym na podstawie otrzymanych wektorów y1i y2

możemy obliczyć x = Qy.

Algorytm z mnożnikami Lagrange’a: Jeśli macierz B jest kolumnowo regularna, to symetryczna macierz M = B^TBjest dodatnio określona. Wtedy kwadrat normy drugiej residuum pierwszego podukładu jest wielomianem drugiego stopnia współrzędnych wektora x, który w każdej warstwie przestrzeni Rⁿma jednoznacznie określone minimum:

f(x) =kd − Bxk²2= x^TMx − 2x^TB^Td + d^Td.

Aby znaleźć minimum funkcji f w zbiorze rozwiązań drugiego podukładu, wystarczy rozwiązać układ równań

"

M C^T

C 0

# "

x y

#

=

"

B^Td e

# .

Po rozwiązaniu tego układu blok y odrzucamy; jego współrzędne są tzw.

mnożnikami Lagrange’a dla postawionego tu zadania minimalizacji z więzami.

Podany wyżej blokowy układ równań możemy rozwiązać w podobny sposób, jak układ w zadaniu 7 na s. 4.26, ale zamiast obliczenia macierzy M i zastosowania do niej metody Choleskiego, lepiej jest dokonać rozkładu ortogonalno-trójkątnego macierzy B.

Opracuj szczegóły obu algorytmów (w szczególności określ wymiary i sposoby reprezentowania wszystkich macierzy otrzymanych w obliczeniach).

Uzasadnij stwierdzenie, że macierz BQ2 przetwarzana w pierwszym algorytmie jest kolumnowo-regularna.

Udowodnij, że drugi algorytm daje rozwiązanie zadania.

5. Opisz, jak zrealizować pierwszy podany na wykładzie algorytm rozwiązywania DLZNK korzystający z rozkładu trójkątno-ortogonalnego, za pomocą odbić Householdera (bez jawnego wyznaczania macierzy Q1).

6. Metoda Newtona z pseudoodwrotnością służy do numerycznego rozwiązywania układów równań nieliniowych, w których liczba równań, m, jest mniejsza niż liczba niewiadomych, n. Układ równań liniowych Jkδ = −fkjest rozwiązywany jako DLZNK, w celu wyznaczenia najkrótszego spełniającego ten układ wektora δ, po czym przyjmuje się xk+1= xk+ δ. W ten sposób, jeśli funkcja f spełnia odpowiednie warunki, powstaje ciąg (xk)k∈Nzbieżny do pewnego

rozwiązania α (z nieskończonego zbioru rozwiązań), położonego w pobliżu punktu startowego x0. Przykład (jedno równanie z dwiema niewiadomymi) na obrazku.

x_k

x_k+1 x

y

xk

xk+1

f₁(x) = 0

x

y z

z = f1(x)

Wykonaj dwa kroki metody Newtona z pseudoodwrotnością dla układu równań

 x²+ y²− z²= 3 x + y + z = 1

przyjmując x0= [0, 0, 0]^T. W rachunkach „ręcznych” układaj i rozwiązuj dualny układ równań normalnych.

7. Udowodnij, że jeśli macierz A ∈ R^m,njest kolumnowo-regularna i jest iloczynem macierzy Q1∈ R^m,ni R1∈ R^n,n, takich że macierz Q1 ma kolumny wzajemnie prostopadłe i o długości 1, zaś macierz R1jest trójkątna górna, to rozkład ten jest jednoznaczny z dokładnością do zwrotów kolumn macierzy Q1i wierszy

macierzy R1.

Z powyższego stwierdzenia wynika wniosek, że wszystkie metody znajdowania rozkładu ortogonalno-trójkątnego (Grama-Schmidta, Householdera, Givensa) dają identyczne wyniki z dokładnością do błędów zaokrągleń i zwrotów kolumn Q1

i wierszy R1 (ale, z wyjątkiem metody Grama-Schmidta, macierzy Q1 nie wyznaczamy w jawnej postaci, jeśli tylko nie musimy!).

Zauważ, że jeśli macierz A jest kwadratowa n × n i ma rząd n − 1 lub n, to jej rozkład na macierz ortogonalną Q i trójkątną górną R też jest jednoznaczny z dokładnością do zwrotów kolumn Q i wierszy R.

8. Sprawdź, że w metodzie sprzężonych gradientów obliczenie r_k+1= r_k− t_kAv_k

jest krokiem ortogonalizacji Grama-Schmidta w sensie iloczynu skalarnego hu, vi = v^Tu, zaś konstrukcja

vk+1= r_k+1+ s_kvk

jest krokiem ortogonalizacji Grama-Schmidta w sensie iloczynu skalarnego hu, vi^A= v^TAu.