Algebra lineare numerica

In questa sezione andiamo ad indagare metodi di scomposizione, iterativi e non. Ci sono molte matrici importanti per questa parte che dovremmo prendere confidenza.

Immagini

![[Lab 2 images]]

Metodo di gauss

Vogliamo cercare un metodo per calcolare soluzioni a sistemi di equazione del genere:

$Ax=b$, classico. Supponiamo che questo sistema abbia una soluzione.

Il nostro obiettivo sarebbe scomporre la matrice $A=LU$ come prodotto di due matrici Lower triangular e Upper triangular.

Algoritmo per matrici Upper e lower triangular

Nel caso io abbia una matrice Lower o Upper triangular, i sistemi nella forma $Lx=b$ sono risolvibili con un algoritmo abbastanza banale che qui non riporto (sostituzioni via via…) Il costo è di $O(n^2)$ e non ti scordare il /2 perché a lei piace (poi ho detto su alla professoressa di Bologna che diceva che $O\left(\frac{n^2}{2}\right)$ è diverso da quanto sopra.)

Scomposizione A = LU

Questo è un algoritmo che necessita di una osservazione abbastanza difficile:

Noto che durante l'applicazione dell’algoritmo di gauss, per creare il pivot su una colonna, è esattamente come se moltiplicassi per una (matrice identità + fattori per la colonna di interesse).

Come in esempio

Si può notare inoltre che che l’inverso di L1 è in una forma molto bella, esattamente la stessa con numeri invertiti di segno!, inoltre la moltiplicazione fra Ln è anch’essa in una forma molto carina! Facilita molto il prodotto. In questo modo mi creo una matrice L, in cui è molto facile invertire.

costo $O(n^3/3)$ per la fattorizzazione LU e poi n2 due volte per risolvere L e U

Upper triangular matrix inversion

Let $\mathbf{U}$ be an $n\times n$ upper triangular matrix. We want to find its inverse $\mathbf{V}={\mathbf{U}}^{-1}$, such that $\mathbf{UV}=\mathbf{I}$, where $\mathbf{I}$ is the $n\times n$ identity matrix. Let $u_{ij}$ be the element in the $i$-th row and $j$-th column of $\mathbf{U}$, and $v_{ij}$ be the element in the $i$-th row and $j$-th column of $\mathbf{V}$.

The algorithm proceeds as follows:

1. Initialize the inverse matrix: Create an $n\times n$ matrix $\mathbf{V}$ and initialize it to the identity matrix $\mathbf{I}$.

2. Iterate through the rows from bottom to top: For each row $i$ from $n$ down to 1:

   • Check for non-zero diagonal element: If $u_{ii}=0$, then the matrix $\mathbf{U}$ is singular and not invertible. Stop the algorithm and indicate this.
   • Normalize the current row: Divide each element in the $i$-th row of $\mathbf{V}$ by $u_{ii}$. This effectively makes the diagonal element of $\mathbf{U}$ in the $i$-th row equal to 1 when considering the product $\mathbf{UV}$.

3. Eliminate off-diagonal elements above the current row: For each row $j$ from $i-1$ down to 1:

   • Calculate the multiplier $m=u_{ji}/u_{ii}$.
   • Subtract $m$ times the $i$-th row of $\mathbf{V}$ from the $j$-th row of $\mathbf{V}$. This will make the element $u_{ji}$ (which is above the diagonal) effectively zero in the context of the inversion process.

4. The resulting matrix $\mathbf{V}$ is the inverse of $\mathbf{U}$.

In essence, the algorithm performs a form of Gaussian elimination, working from the bottom right of the matrix upwards and towards the left. Because of the upper triangular structure, we only need to eliminate the elements above the diagonal.

Example:

Let's consider the upper triangular matrix: $\mathbf{U}=\left(\begin{array}{ccc}2& 1& 3\\ 0& -1& 2\\ 0& 0& 4\end{array}\right)$

1. Initialize $\mathbf{V}=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right)$.

2. Row 3 ($i=3$):

   • $u_{33}=4\ne 0$.
   • Normalize row 3 of $\mathbf{V}$: $\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0/4& 0/4& 1/4\end{array}\right)=\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1/4\end{array}\right)$.

3. Row 2 ($i=2$):

   • $u_{22}=-1\ne 0$.
   • Normalize row 2 of $\mathbf{V}$: $\left(\begin{array}{ccc}1& 0& 0\\ 0/-1& 1/-1& 0/-1\\ 0& 0& 1/4\end{array}\right)=\left(\begin{array}{ccc}1& 0& 0\\ 0& -1& 0\\ 0& 0& 1/4\end{array}\right)$.
   • Eliminate above: $j=1$. $m=u_{12}/u_{22}=1/(-1)=-1$.
     • Row 1 of $\mathbf{V}$ becomes: $\left(\begin{array}{ccc}1-(-1)\times 0& 0-(-1)\times (-1)& 0-(-1)\times 0\end{array}\right)=\left(\begin{array}{ccc}1& -1& 0\end{array}\right)$.
     • $\mathbf{V}=\left(\begin{array}{ccc}1& -1& 0\\ 0& -1& 0\\ 0& 0& 1/4\end{array}\right)$.

4. Row 1 ($i=1$):

   • $u_{11}=2\ne 0$.
   • Normalize row 1 of $\mathbf{V}$: $\left(\begin{array}{ccc}1/2& -1/2& 0/2\\ 0& -1& 0\\ 0& 0& 1/4\end{array}\right)=\left(\begin{array}{ccc}1/2& -1/2& 0\\ 0& -1& 0\\ 0& 0& 1/4\end{array}\right)$.
   • Eliminate above: (no rows above row 1).

Therefore, the inverse of $\mathbf{U}$ is: ${\mathbf{U}}^{-1}=\left(\begin{array}{ccc}1/2& -1/2& 0\\ 0& -1& 0\\ 0& 0& 1/4\end{array}\right)$

You can verify this by multiplying $\mathbf{U}$ and ${\mathbf{U}}^{-1}$ to see if you get the identity matrix.

This algorithm efficiently computes the inverse by taking advantage of the zero entries in the lower triangle, reducing the number of operations required.

Scomposizione con permutazione

Ogni tanto potrebbe capitare che ci sia il bisogno di permutare, altrimenti non ho nessuna soluzione, da notare che questo è un passo all'interno dell’algoritmo di gauss.

Inoltre scegliere il maggiore dovrebbe aiutare ad eliminare mini-errorini dovuti all’imprecisione della somma floating point.

• Algoritmo

  

Forma finale scomposizioni

Se utilizziamo una sorta di pivoting parziale ossia per ridurre gli errori scegliamo il massimo della colonna (come del resto fa l’algoritmo di sopra) allora non abbiamo bisogno di permutare colonne), se vogliamo invece un pivoting totale allora dobbiamo permutare le colonne e la matrice diventa qualcosa del tipo

$A=QLUP$ con Q le permutazioni su riga, e P quelle su colonna (quelle su colonna non è che ci importa molto, perché basta moltiplicare per $P^T$, ossia la sua inversa per x e si risolvono problemi di questo genere… dato che alla fine $P{P}^T=I$

Metodo di Cholesky

Se prendo una matrice A simmetrica definita positiva (ricorda analisi per questo), allora con questo metodo riesco a riscriverla nella forma $A=L{L}^T$, e so che la matrice L triangolare inferiore non è singolare perché ho tutti gli autovalori positivi.

Questo metodo costa $O(n^3/6)$ leggermente più veloce della fattorizzazione di Gauss.

Una volta fatto la scomposizione poi va subito

$Ly=b,L^{T}x=y$ e si risolve, sul come funziona l'algoritmo, per ora ci importa solamente che è più veloce.

Se non è simmetrico o non è definita positiva il programma crasha.

Metodo di Schur

Questo è un metodo nuovo che aggiungo in seguito. Una risorsa buona per questo metodo lo trovi qui. È utile quando vogliamo fare una cosa di questo tipo

$$A=QB{Q}^{-1}$$

Dove $Q$ è una matrice di cambio di base. Questa cosa non si può fare sempre quando vogliamo che $B$ sia diagonale, però la decomposizione di Schur ci dà qualche nota in più riguardo a questo, usando però matrici upper triangular

Enunciato di Schur

Per qualunque matrice a valori complessi $A$ esiste una base ortonormale in $\mathbb{C}$ e una matrice upper triangular $T$ tale per cui vale

$$A=BT{B}^{-1}$$

NOTE: se $A$ è una matrice normale, allora $T$ è diagonale. Normale quando vale $A{A}^{\ast }=A^{\ast }A$. Con $A^{\ast }$ la trasposta coniugata della matrice $A$.

Idee generali per la dimostrazione di Schur

In pratica per un teorema dell'algebra (quello fondamentale) si ha che $T$, dato che è quadrata, ha almeno un autovettore, allora prendi questo autovettore, usi il teorema di estensione per creare una base, e scomponi con questo. Questo ti creerà una matrice a blocchi grossi $2\times 2$ in cui quello in basso a destra è ancora quadrato con suo sottospazio. Continui così finché non finisce tutto e resta una matrice upper triangular. Nella reference ciò è descritto molto bene da un punto di vista intuitivo. La cosa carina con le $T$ è che per esempio è spesso più facile elevarli alla potenza, non quanto le diagonali, ma più semplice.

Matrici semi-definite o definite positive

Se prendo un qualunque vettore e vale

$$x^{T}Ax\ge 0$$

Oppure avendo autovalori sempre positivi sono maggiore di 0, si può vedere che è equivalente. Certi autori richiedono anche che la matrice $A$ sia simmetrica.

Guarda Massimi minimi multi-variabile nella sezione forme quadratiche per questo.

Metodi iterativi

Sono buoni perché l'errore di questi è molto più alto, ma spesso molto più veloce, mentre per i metodi diretti è esattamente il contrario (errore basso, alta complessità).

In modo molto generale, possiamo definire un metodo iterativo come una funzione $G$ applicata all'input stesso: $x_n=G(x_{n-1})$

Note introduttive convergenza e iterazione

Questa sotto è la definizione di convergenza ad $\alpha$ con ordine $p\ge 1$

Spesso la convergenza è fissata da una condizione di arresto (che può essere dipendente anche dal numero di iterazioni. Possiamo definire C = fattore di convergenza quando $p=1$, perché più è piccolo più converge in fretta. C'è anche una relazione simile fra i p, più è grande, in questo caso, più velocemente converge.

Metodi iterativi stazionari

Idee Generali

Questi sono nella forma $x_{k+1}=H{x}_k+d_k$ Poi scriverlo come una differenza di matrici in particolare vorre $A=M-N$ , dove $M$ sia non singolare (e possibilmente facilmente invertibile). Se ciò è possibile allora

$Ax=b⟹Mx-Nx=b⟹x=M^{-1}Nx+M^{-1}b⟹x=Tx+c$

Una volta avuto questa matrice $T$ possiamo riutilizzarla per ogni iterazione di successione!

Avremo che $x_k=T{x}_{k-1}+c$, e si ha che $x\ast =\lim x_k$ deve essere che

$$x\ast =Tx\ast +c$$

E vorremmo che questo metodo converga per ogni input quindi potrei mettere un x_0 a casissimo

Una cosa non espressa durante la lezione è che $N=M-A$ quindi

$$M^{-1}N=M^{-1}(M-A)=I-M^{-1}A$$

Che ci permette di scrivere la stessa cosa sente fare utilizzo della matrice $N$ e secondo me è più clean. Questa cosa è presente nel libro.

Condizione necessaria e sufficiente di convergenza

Teorema sulla condizione di convergenza:

Questo risultato si può connettere col raggio spettrale e si può concludere che converge sse

$$\rho (T)<1$$

Metodo di Jacobi

Sia $D$ la parte diagonale di $A$, allora l’iterazione in questa forma converge e si chiama metodo di jacobi

$$x_{k+1}=(I-D^{-1}A)x_k+D^{-1}b$$

Nota: attento che la prof non ha insegnato in questo modo! Quindi probabilmente se te la chiede diglielo nella forma che le piace,

$M=D,N=E+F$

Metodo Gauss-Seidel

Si ricorda che $A=D-E-F$ in modo molto strano di scriverlo…

Si fa uso della matrice Lower con anche la diagonale, per il resto è esattamente uguale a tutto il resto per i metodi iterativi., quindi

$$M=D-E,N=F$$

Condizioni di convergenza per Jacobi e Gauss-Seidel

Per comprendere questa parte è importante il concetto di matrice con diagonale dominante. In parole umane deve essere che ogni elemento sulla matrice deve essere maggiore della somma di tutto quanto non stia sulla diagonale maggiore.

• Slide condizioni

  

Nota questi sono in pratica da imparare a memoria, ma non me ne ricordo!

Gradiente coniugato

Questo è un metodo di arresto molto misterioso, direi di cacharla 🤑.

Metodi di cui non ha parlato

Metodo SOR

Metodi di Rilassamento

Metodi di Krylov