1. naloga

(a) Ker je rang razširjeme matrike sistema \([A | \mathbf{b}]\) enak \(3\), rang \(A\) pa enak \(2\), sistem \(A\mathbf{x} = \mathbf{b}\) nima rešitev. Pravokotno projekcijo \(\mathbf{b}\) na \(C(A)\) bomo dobili tako, da poiščemo rešitev \(\mathbf{x}\) v smislu najmanjših kvadratov za \(A \mathbf{x} = \mathbf{b}\), potem je \(\mathbf{b}_1 = A\mathbf{x}\). Rešitev \(\mathbf{x}\) po metodi najmanjših kvadratov najlažje poiščemo kot rešitev normalnega sistema \(A^\mathsf{T} A \mathbf{x} = A^\mathsf{T} \mathbf{b}\). V našem primeru je \[A^\mathsf{T} A = \begin{bmatrix} 4 & -4 & 0 \\ -4 & 4 & 0 \\ 0 & 0 & 16 \end{bmatrix} \ , \ A^\mathsf{T} \mathbf{b} = \begin{bmatrix} -4 \\ 4 \\ 16 \end{bmatrix}.\] Gaussova eliminacija na \([A^\mathsf{T} A | A^\mathsf{T} \mathbf{b}]\) pa da \[\left[\begin{array}{ccc|c} 1 & -1 & 0 & -1 \\ 0 & 0 & 1 & 1 \\ 0 & 0 & 0 & 0 \end{array}\right].\] Sledi, da je splošna rešitev \(\mathbf{x} = [x_2 - 1, x_2, 1]^\mathsf{T}\), kjer je \(x_2\) poljubno realno število. Pravokotna projekcija \(\mathbf{b}\) na \(C(A)\) je neodvisna od \(x_2\): \[\mathbf{b}_1 = A\mathbf{x} = \begin{bmatrix} 3 \\ -3 \\ 1 \\ -1 \end{bmatrix}.\] 

Sistem \(A'\mathbf{x} = \mathbf{b}'\) prav tako nima rešitev, rešitev pripadajočega normalnega sistema pa je \(\mathbf{x} = [-1,1]^\mathsf{T}\). Pripadajoča pravokotna projekcija je enaka kot prej; \(\mathbf{b}_1' = \mathbf{b}_1\), kar bi lahko opazili hitreje, saj sta stolpčna prostora matrik \(A\) in \(A'\) enaka.

(b) Poiskali bomo t.i. ekonomični singularni razcep \(A = USV^\mathsf{T}\), v katerem je \(S\) kvadratna diagonalna matrika s singularnimi vrednostmi na diagonali, ena od ortogonalnih matrik \(U\) in \(V\) pa je pravokotna. Iz \(A = USV^\mathsf{T}\) sledi \(A^\mathsf{T} A = VS^2 V^\mathsf{T}\) (in tudi \(AA^\mathsf{T}  = US^2 U^\mathsf{T}\)). (Zakaj?) To pomeni, da so kvadrati singularnih vrednosti matrike \(A\) ravno lastne vrednosti matrike \(A^\mathsf{T}A\), tj. \(\lambda_i = \sigma_i^2\) za lastno vrednost \(\lambda_i\) matrike \(A^\mathsf{T} A\) in singularno vrednost \(\sigma_i\) matrike \(A\), stolpci matrike \(V\) pa so ravno normirani lastni vektorji matrike \(A^\mathsf{T} A\). Poiščimo najprej lastne vrednosti matrike \(A^\mathsf{T} A\), te so ničle karakterističnega polinoma \[\det(A^\mathsf{T}A-\lambda I) = \begin{vmatrix} 4-\lambda & -4 & 0 \\ -4 & 4-\lambda & 0 \\ 0 & 0 & 16-\lambda \end{vmatrix} = -(16-\lambda)(8-\lambda)\lambda.\] Lastne vrednosti \(A^\mathsf{T}A\) so torej \(\lambda_1=16\), \(\lambda_2 = 8\) in \(\lambda_3 = 0\). Matrika \(S\) je torej \[S = \begin{bmatrix}\sqrt{\lambda_1} & 0 & 0 \\ 0 & \sqrt{\lambda_2} & 0 \\ 0 & 0 & \sqrt{\lambda_3} \end{bmatrix} = \begin{bmatrix}4 & 0 & 0 \\ 0 & 2\sqrt{2} & 0 \\ 0 & 0 & 0 \end{bmatrix}.\] Pripadajoče lastne vektorje dobimo kot (neničelne) rešitve homogenih sistemov \((A^\mathsf{T}A - \lambda_i I)\mathbf{v} = \mathbf{0}\), ti so \[\mathbf{v}_1 = \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix}, \mathbf{v}_2 = \begin{bmatrix} 1 \\ -1 \\ 0 \end{bmatrix}, \mathbf{v}_3 = \begin{bmatrix} 1 \\ 1 \\ 0 \end{bmatrix}.\] Prvi je že normiran, druga dva pa ne, zato ju normiramo, preden jih vpišemo v matriko \(V\): \[V = \big[\mathbf{v}_1, \frac{1}{\sqrt{2}}\mathbf{v}_2, \frac{1}{\sqrt{2}}\mathbf{v}_3\big] = \frac{1}{\sqrt{2}}\begin{bmatrix} 0 & 1 & 1 \\ 0 & -1 & 1 \\ \sqrt{2} & 0 & 0\end{bmatrix}.\] Poiskati moramo še matriko \(U\). Iz \(A = USV^\mathsf{T}\) sledi \(AV = US\). Iščemo katerokoli ortogonalno matriko \(U\), da velja \[ AV = \begin{bmatrix} 2 & -\sqrt{2} & 0 \\ -2 & \sqrt{2} & 0 \\ 2 & \sqrt{2} & 0 \\ -2 & -\sqrt{2} & 0 \end{bmatrix} = U \begin{bmatrix}4 & 0 & 0 \\ 0 & 2\sqrt{2} & 0 \\ 0 & 0 & 0 \end{bmatrix} = US. \] Prva dva stolpca \(U\) sta torej povsem določena, tretjega pa lahko izberemo skoraj poljubno - matrika \(U\) mora biti ortogonalna. Dobimo (npr.) \[ U = \frac{1}{2}\begin{bmatrix} 1 & -1 & \sqrt{2} \\ -1 & 1 & \sqrt{2} \\ 1 & 1 & 0 \\ -1 & -1 & 0 \end{bmatrix}.\] Ničelne singularne vrednosti lahko zanemarimo, saj nič ne prispevajo k matriki \(A\). Če jih, dobimo ti. kompaktni oz. tanki singularni razcep matrike \(A\), \(A = U' S' V'^\mathsf{T}\), kjer je \[U' = \frac{1}{2}\begin{bmatrix} 1 & -1\\ -1 & 1 \\ 1 & 1 \\ -1 & -1 \end{bmatrix}, S' = \begin{bmatrix}4 & 0 \\ 0 & 2\sqrt{2} \end{bmatrix}, V' = \frac{1}{\sqrt{2}}\begin{bmatrix} 0 & 1 \\ 0 & -1 \\ \sqrt{2} & 0 \end{bmatrix}.\] (Iz razcepa torej izpustimo vse ničelne singularne vrednosti ter pripadajoče leve in desne singularne vektorje.)

(c) Iz singularnega razcepa \(A = USV^\mathsf{T}\) Moore-Penroseov posplošeni inverz izračunamo po formuli \(A^+ = V S^+ U^\mathsf{T}\), kjer so na diagonali matrike \(S^+\) vrednosti \(1/\sigma_i\), če je \(\sigma_i > 0\), in ničle za vsako ničelno singularno vrednost. V resnici lahko vedno uporabimo kar tanki SVD za \(A\), tedaj je \(S'^+ = S'^{-1}\) in \[A^+ = V' S'^{-1} U'^\mathsf{T} = \frac{1}{8} \begin{bmatrix} -1 & 1 & 1 & -1 \\ 1 & -1 & -1 & 1 \\ 1 & -1 & 1 & -1 \end{bmatrix}.\] Od tod dobimo \[A^+ \mathbf{b} = \frac{1}{2}\begin{bmatrix} -1 \\ 1 \\ 2 \end{bmatrix}.\] To pa je točno tista rešitev \(\mathbf{x}\) iz (b) dela naloge, ki ima najmanjšo normo. (Preveri to!) 

Ker ima matrika \(A'\) poln stolpčni rang, lahko njen Moore-Penroseov posplošeni inverz izračunamo po formuli \(A'^+ = (A'^\mathsf{T}A')^{-1} A'^\mathsf{T}\). Dobimo \[A'^+ = \frac{1}{8} \begin{bmatrix} -2 & 2 & 2 & -2 \\ 1 & -1 & -1 & 1 \end{bmatrix}.\] Velja seveda \(A'^+\mathbf{b} = [-1,1]^\mathsf{T}\), kar je ravno rešitev sistema \(A'\mathbf{x} = \mathbf{b}\) v smislu linearne metode najmanjših kvadratov, ki smo jo poiskali v (a) delu naloge.

(d) Najprej shranimo matriko \(A\) in vektor \(\mathbf{b}\):

octave:1> A = [-1 1 2; 1 -1 -2; 1 -1 2; -1 1 -2]
octave:2> b = [3; -3; 2; 0]

Klic

octave:3> [U, S, V] = svd(A)

poišče polni singularni razcep matrike \(A\). Ekonomični SVD dobimo z:

octave:4> [U, S, V] = svd(A, "econ")

Če smo pozorni, opazimo, da je tretja singularna vrednost zelo majhna (reda velikosti strojne natančnosti), ni pa \(0\). Če bi naivno zaupali formuli, bi za Moore-Penroseov posplošeni inverz dobili

octave:5> Aplus = V*inv(S)*U'

kar je seveda daleč od našega rezultata. Singularne vrednosti, ki so le skoraj \(0\) zaradi zaokrožitvenih napak pri računih v plavajoči vejici, moramo interpretirati kot \(0\). Klic 

octave:6> Aplus = pinv(A)

to upošteva in vrne enak rezultat kot smo ga dobili zgoraj.

Zadnja sprememba: nedelja, 10. marec 2024, 20.55