Rešitve 1. naloge, 3. teden

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.