Algebra lineare numerica^{1}

A. Sommariva^{2}

Keywords: ...

1. Condizionamento di matrici e sistemi

Definizione. Dato uno spazio vettoriale V, la funzione ∥ · ∥ : V → ℝ è una norma vettoriale se

1. ∥x∥ ≥ 0 per ogni x ∈ V e ∥x∥ = 0 se e solo se x = 0_{V},

2. ∥αx∥ = |α|∥x∥ per ogni x ∈ V e scalare α ∈ ℝ,

3. vale la cosidetta disuguaglianza triangolare

∥x + y∥ ≥ ∥x∥ + ∥y∥

per ogni x, y ∈ V.

Esempio. Se V = ℂ^{n} e x = (x_{1}, x_{2}, …, x_{n}) ∈ ℂ^{n} alcuni esempi di norme sono

• per p ∈ [1, +∞), la norma p è definita come

∥x∥_{p} = (∑_{j=1}^{n} |x_{j}|^{p})^{1/p};

• la norma del massimo è definita come

∥x∥_{p} = max_{j=1,...,n} |x_{j}|;

Nota. 1.1. Lo stesso tipo di norme si definisce nel caso particolare i vettori appartengano a ℝ^{n}.

Nota. 1.2. Si osservi che

• per p = 1, si definisce la importante norma 1

∥x∥_{1} = ∑_{j=1}^{n} |x_{j}|;

• per p = 2, si definisce la importante norma 2, detta anche euclidea

∥x∥_{2} = (∑_{j=1}^{n} |x_{j}|^{2})^{1/2} .

Esempio. Si consideri il vettore v = (4, −3) ∈ ℝ^{2}. Allora • ∥v∥_{1} = |4| + | − 3| = 7,

• ∥v∥_{2} = \sqrt{|4|^{2} + | − 3|^{2}} = 5,

• ∥v∥_{∞} = max |4|, | − 3| = 4. Esempio. Si considerino i vettori u = (4, −3), v = (4.001, −2.999) ∈ ℝ^{2}.

I due vettori sono numericamente vicini. Per verificarlo calcoliamo la loro distanza in termine di norme ovvero

dist(u, v) := ∥u − v∥_{p}, p = [1, ∞].

Allora, visto che u − v = [−0.001 − 0.001]

• ∥u − v∥_{1} = 0.002,

• ∥u − v∥_{2} = 0.0014,

• ∥u − v∥_{∞} = max |4|, | − 3| = 0.001, risultati che sottolineano che u e v sono vicini rispetto alle distanze definite dalle norme. Definizione. Sia ∥·∥ una norma vettoriale. Si definisce norma vettoriale indotta di una matrice in ℂ^{n×n} la funzione che mappa ogni matrice A ∈ ℝ^{n×n} nel valore

∥A∥ = sup_{x≠0} \frac{∥Ax∥}{∥x∥}.

Nota. 1.3. Lo stesso tipo di norme si definisce nel caso particolare le matrici siano in ℝ^{n×n}.

Alcuni esempi di norme matriciali risultano:

• per p = 1, partendo dalla norma vettoriale ∥ · ∥_{1}, si dimostra che la norma 1 matriciale risulta uguale a

∥A∥_{1} = max_{1≤j≤n} ∑_{i=1}^{n} |a_{i,j}|,

dove A = (a_{i,j}), • per p = 2, partendo dalla norma vettoriale ∥ · ∥_{2}, si dimostra che la norma 2

matriciale risulta uguale a

∥x∥_{2} = \sqrt{σ(A^{∗}A)}

dove A = (a_{i,j}), σ(B) è il raggio spettrale di B ovvero il massimo modulo di un autovalore di B, e se A ∈ ℂ^{n×n} allora A^{∗} è la trasposta coniugata di A mentre se A ∈ ℝ^{n×n} allora A^{∗} è la trasposta di A;

• per p = ∞, partendo dalla norma vettoriale ∥ · ∥_{∞}, si dimostra che la norma ∞ matriciale risulta uguale a

∥A∥_{∞} = max_{1≤i≤n} ∑_{j=1}^{n} |a_{i,j}|,

dove A = (a_{i,j}).

Esempio. Si consideri la matrice A ∈ ℝ^{2×2}

A = [ \begin{matrix}1 & -2 \\ -3 & 4 \end{matrix} ]

Allora

• ∥A∥_{1} = max(|1| + | − 3|, | − 2| + |4|) = 6;

• osservato che A^{∗} = A^{T} abbiamo

A^{T} A = [ \begin{matrix}1 & -3 \\ -2 & 4 \end{matrix} ] [ \begin{matrix}1 & -2 \\ -3 & 4 \end{matrix} ] = [ \begin{matrix}10 & -14 \\ -14 & 20 \end{matrix} ]

i cui autovalori valgono rispettivamente

λ_{1} ≈ 0, 1339312526814949, λ_{2} ≈ 29.86606874731851

e di conseguenza

∥A∥_{2} = \sqrt{max(|λ_{1}|, |λ_{2}|)} ≈ 5.464985704219044;

• ∥A∥_{∞} = max(|1| + | − 2|, | − 3| + |4|) = 7.

Esempio. Si considerino le matrici

A = [ \begin{matrix}1.001 & -2.999 \\ -2.002 & 4.005 \end{matrix} ] B = [ \begin{matrix}1 & -3 \\ -2 & 4 \end{matrix} ]

Le due matrici sono numericamente vicine. Per verificarlo calcoliamo la loro distanza in termine di norme ovvero

dist(A, B) := ∥A − B∥_{p}, p = [1, ∞].

Allora, visto che

A − B = [ \begin{matrix}0.0010 & 0.0010 \\ -0.0020 & 0.0050 \end{matrix} ]

ricaviamo • ∥A − B∥_{1} = 0.006,

• ∥A − B∥_{2} = 0.005415654779058332,

• ∥A − B∥_{∞} = max (|4|, | − 3|) = 0.007,

risultati che sottolineano che A e B sono matrici vicine rispetto alle distanze definite dalle norme.

Teorema 1.4. Per le norme matriciali indotte, valgono le seguenti disuguaglianze

∥Ax∥ ≤ ∥A∥∥x∥,  (1)

e

∥AB∥ ≤ ∥A∥∥B∥.  (2)

Dimostrazione. Mostriamo la prima disuguaglianza. E’ banalmente verificata per x = 0. Altrimenti, essendo

∥A∥ = sup_{y≠0} \frac{∥Ay∥}{∥y∥} ≥ \frac{∥Ax∥}{∥x∥}

abbiamo che per x ≠ 0,

∥Ax∥ ≤ ∥A∥∥x∥.

Per quanto riguarda la seguente disuguaglianza, essendo per il primo punto, se y ≠ 0

∥ABy∥ = ∥A(By)∥ ≤ ∥A∥∥By∥ ≤ ∥A∥∥B∥∥y∥

da cui

\frac{∥ABy∥}{∥y∥} ≤ ∥A∥∥B∥

e quindi

∥AB∥ = sup_{y≠0} \frac{∥ABy∥}{∥y∥} ≤ ∥A∥∥B∥.

Definizione. Sia A ∈ ℝ^{n×n} una matrice invertibile ossia tale che det(A) ≠ 0. Si dice indice di condizionamento della matrice A, relativamente alla norma indotta ∥ · ∥,

k(A) := ∥A∥∥A^{−1}∥.

Nota. 1.5. Si osservi che dalla definizione, se I è la matrice identica di ℝ^{n×n}

∥I∥ = sup_{x≠0} \frac{∥Ix∥}{∥x∥} = sup_{x≠0} \frac{∥x∥}{∥x∥} = 1

e quindi

1 = ∥I∥ = ∥AA^{−1}∥ ≤ ∥A∥∥A^{−1}∥ = k(A)

ovvero k(A) ≥ 1.

Vale il seguente risultato

Teorema 1.6. Qualsiasi sia la norma matriciale indotta, e A invertibile,

1 ≤ \frac{|λmax|}{|λmin|} ≤ k(A)

dove λmin e λmax sono rispettivamente gli autovalori di A di minimo e massimo modulo.

Se A è simmetrica

k_{2}(A) = \frac{|λmax|}{|λmin|}.

Esempio. La matrice H = (h_{i,j})_{i,j=1,...,n} con h_{i,j} = \frac{1}{i+j−1} è nota come matrice di Hilbert di ordine n. Per queste matrici gli indici di condizionamento crescono molto rapidamente, come si può vedere dalla tabella.

\begin{matrix}n & c_{n} & n & c_{n} & n & c_{n} \\ 1 & 1.00e + 00 & 5 & 4.77e + 05 & 9 & 4.93e + 11 \\ 2 & 1.93e + 01 & 6 & 1.50e + 07 & 10 & 1.60e + 13 \\ 3 & 5.24e + 02 & 7 & 4.75e + 08 & 11 & 5.22e + 14 \\ 4 & 1.55e + 04 & 8 & 1.53e + 10 & 12 & 1.62e + 16 \end{matrix}

Tabella 1: Indice di condizionamento c_{n} = cond_{2}(H_{n}) della matrici di Hilbert H_{n} di ordine n

Vale la seguente risposta della soluzione di un sistema lineare agli errori del vettore termine noto.

Teorema 1.7. Se

1. A ∈ ℝ^{n×n} è una matrice invertibile,

2. b ∈ ℝ^{n}, b ≠ 0,

3. Ax = b,

4. A(x + δx) = b + δb,

allora vale la stima

\frac{∥δx∥}{∥x∥} ≤ k(A)\frac{∥δb∥}{∥b∥}.

Dimostrazione. Da A(x + δx) = b + δb e Ax = b abbiamo che

Aδx = A(x + δx) − Ax = (b + δb) − b

ovvero Aδx = δb e quindi δx = A^{−1}δb. Ne consegue che

∥δx∥ = ∥A^{−1}δb∥ ≤ ∥A^{−1}∥∥δb∥.  (3)

D’altra parte, essendo Ax = b

∥b∥ = ∥Ax∥ ≤ ∥A∥∥x∥

e quindi, essendo x ≠ 0 in quanto b ≠ 0, ricaviamo

\frac{1}{∥x∥} ≤ \frac{∥A∥}{∥b∥}.  (4)

Da (3), (4), ricaviamo

\frac{∥δx∥}{∥x∥} ≤ ∥A^{−1}∥∥A∥\frac{∥δb∥}{∥b∥} = k(A)\frac{∥δb∥}{∥b∥}.

Commento. Questo teorema mostra che si fa un errore relativo sul termine noto pari a \frac{∥δb∥}{∥b∥} allora si compie un errore relativo sulla soluzione \frac{∥δx∥}{∥x∥} tale che

\frac{∥δx∥}{∥x∥} ≤ k(A)\frac{∥δb∥}{∥b∥}.

Quindi, se k(A) ≫ 1, può accadere che la soluzione perturbata x + δx sia molto distante da x e di conseguenza la soluzione del sistema lineare Ax = b molto suscettibile a piccoli errori sul termine noto. In altri termini più grande è k(A) più i sistemi lineari Ax = b sono difficili da trattare nel senso che piccoli errori sui dati b possono portare a grossi errori sulla soluzione x.

Esempio. Siano

A = [ \begin{matrix}7 & 10 \\ 5 & 7 \end{matrix} ] , b = [ \begin{matrix}1 \\ 0.7 \end{matrix} ] , b + δb = [ \begin{matrix}1.01 \\ 0.69 \end{matrix} ] .

Il sistema Ax = b ha soluzione

x = [ \begin{matrix}0 \\ 0.1 \end{matrix} ] .

mentre Ax = b + δb ha soluzione

x + δx = [ \begin{matrix}-0.17 \\ 0.22 \end{matrix} ] .

Di conseguenza seppure \frac{∥δb∥_{∞}}{∥b∥_{∞}} = 0.01 abbiamo \frac{∥δx∥_{∞}}{∥x∥_{∞}} = 1.7. Questo significa che a un piccolo errore relativo sui dati abbiamo riscontrato un eccessivo errore relativo sulla soluzione.

Osserviamo che in questo caso k(A) = ∥A∥_{∞}∥A^{−1}∥_{∞} = 289 e quindi non troppo vicino a 1.

Mostriamo ora la risposta della soluzione di un sistema lineare agli errori sulla matrice.

Teorema 1.8. Se 1. A ∈ ℝ^{n×n} è una matrice invertibile,

2. b ∈ ℝ^{n}, b ≠ 0,

3. Ax = b,

4. (A + δA)(x + δx) = b, con det(A + δA) ≠ 0, allora vale la stima

\frac{∥δx∥}{∥x + δx∥} ≤ k(A)\frac{∥δA∥}{∥A∥}.

Dimostrazione. Essendo (A + δA)(x + δx) = b e Ax = b, necessariamente

0 = (A + δA)(x + δx) − b = Ax + Aδx + δAx + δAδx − b
= Aδx + δAx + δAδx = Aδx + δA(x + δx)

da cui

Aδx = −δA(x + δx)

ovvero

δx = −A^{−1}δA(x + δx)

e di conseguenza

∥δx∥ = ∥ − A^{−1}δA(x + δx)∥ ≤ ∥A^{−1}∥∥δA∥∥x + δx∥.

Quindi, essendo ∥A∥ ≠ 0,

\frac{∥δx∥}{∥x + δx∥} ≤ ∥A^{−1}∥∥δA∥ = \frac{∥A^{−1}∥∥A∥}{∥A∥}∥δA∥

da cui

\frac{∥δx∥}{∥x + δx∥} ≤ k(A)\frac{∥δA∥}{∥A∥}.

Vogliamo ora studiare il caso in cui tanto la matrice quanto il termine noto siano soggetti ad errori.

Per prima cosa introduciamo il seguente importante teorema. Teorema 1.9. Se

1. ∥ · ∥ è una norma matriciale indotta da una vettoriale,

2. ∥A∥ < 1, allora ”I − A” e ”I + A” sono invertibili e si ha

∥(I − A)^{−1}∥ ≤ \frac{1}{1 − ∥A∥},

∥(I + A)^{−1}∥ ≤ \frac{1}{1 − ∥A∥}.

Finalmente,

Teorema 1.10. Se

1. A ∈ ℝ^{n×n} è una matrice invertibile,

2. b ∈ ℝ^{n}, b ≠ 0,

3. Ax = b,

4. (A + δA)(x + δx) = b + δx, con det(A + δA) ≠ 0,

5. k(A)∥δA∥ < ∥A∥,

allora vale la stima

\frac{∥δx∥}{∥x∥} ≤ \frac{k(A)}{1 − k(A)∥δA∥/∥A∥} (\frac{∥δA∥}{∥A∥} + \frac{∥δb∥}{∥b∥}) .

Commento. Il teorema stima l’errore relativo \frac{∥δx∥}{∥x∥} relativamente alle soluzioni, in funzione degli errori relativi sulla perturbazione della matrice \frac{∥δA∥}{∥A∥} e sul termine noto \frac{∥δb∥}{∥b∥}. Il valore

\frac{k(A)}{1 − k(A)∥δA∥/∥A∥}

è la costante di amplificazione degli errori.

Dimostrazione. Essendo (A + δA)(x + δx) = b + δx e Ax = b deduciamo

b + δb = (A + δA)(x + δx) = Ax + Aδx + δAx + δAδx
= b + Aδx + δAx + δAδx

da cui sottraendo b ad ambo i membri

δb = Aδx + δAx + δAδx = δAx + (A + δA)δx

ovvero

(A + δA)δx = δb − δAx.  (5)

Posto B = A^{−1}δA, ricaviamo

A + δA = A + AA^{−1}δA = A(I + A^{−1}δA) = A(I + B)

e da k(A)∥δA∥ < ∥A∥

∥B∥ = ∥A^{−1}δA∥ ≤ ∥A^{−1}∥∥δA∥
= \frac{∥A^{−1}∥∥A∥}{∥A∥}∥δA∥ = k(A)\frac{∥δA∥}{∥A∥}
≤ \frac{∥A∥}{∥A∥} = 1  (6)

ovvero ∥B∥ < 1.

Osserviamo che la matrice A + δA è invertibile in quanto lo è A, lo è I + B per il teorema 1.9 e da

A + δA = A(I + A^{−1}δA) = A(I + B)

lo è pure A + δA perchè il prodotto di matrici invertibili è invertibile (ricordare che det(MN) = det(M) det(N), e quindi se det(M), det(N) ≠ 0 allora det(MN) ≠ 0).

Dal teorema 5, moltiplicando ambo i membri per (A+δA)^{−1}, visto che (MN)^{−1} = N^{−1}M^{−1} ricaviamo

δx = (A + δA)^{−1}(δb − δAx) = (A(I + B))^{−1}(δb − δAx)
= (I + B)^{−1}A^{−1}(δb − δAx)

(7)

da cui, essendo nuovamente per il teorema 5∥(I+B)^{−1}∥ ≤ 1/(1−∥B∥), B = A^{−1}δA

∥δx∥ = ∥(I + B)^{−1}A^{−1}(δb − δAx)∥
≤ ∥(I + B)^{−1}∥∥A^{−1}∥∥δb − δAx∥
≤ \frac{1}{1 − ∥B∥}∥A^{−1}∥∥δb − δAx∥
= \frac{∥A^{−1}∥}{1 − ∥A^{−1}δA∥}∥δb − δAx∥
≤ \frac{∥A^{−1}∥}{1 − ∥A^{−1}∥∥δA∥}(∥δb∥ + ∥δA∥∥x∥)  (8)

e quindi riassumendo

∥δx∥ ≤ \frac{∥A^{−1}∥}{1 − ∥A^{−1}∥∥δA∥}(∥δb∥ + ∥δA∥∥x∥).  (9)

Essendo b = Ax abbiamo ∥b∥ = ∥Ax∥ ≤ ∥A∥∥x∥ da cui

\frac{1}{∥A∥∥x∥} ≤ \frac{1}{∥b∥}  (10)

e moltiplicando membro a membro i risultati di (9), (10),

\frac{∥δx∥}{∥x∥} ≤ \frac{∥A^{−1}∥}{1 − ∥A^{−1}∥∥δA∥} (\frac{∥δb∥}{∥x∥} + ∥δA∥∥x∥\frac{1}{∥x∥})
≤ \frac{∥A^{−1}∥∥A∥}{1 − ∥A^{−1}∥∥δA∥} \frac{∥1∥}{∥A∥} (\frac{∥δb∥}{∥x∥} + ∥δA∥)
≤ \frac{k(A)}{1 − ∥A^{−1}∥∥A∥∥\frac{∥δA∥}{∥A∥}} (\frac{∥δb∥}{∥A∥∥x∥} + \frac{∥δA∥}{∥A∥})

≤ \frac{k(A)}{1 − k(A)∥\frac{∥δA∥}{∥A∥}} (\frac{∥δb∥}{∥b∥} + \frac{∥δA∥}{∥A∥}) .  (11)

Esempio. Il sistema Ax = b, dove

A = [ \begin{matrix}1 & 2 \\ 0.499 & 1.001 \end{matrix} ] , b = [ \begin{matrix}3 \\ 1.5 \end{matrix} ] ,

ha soluzione x = (1, 1)^{T}.

Se consideriamo

A + δA = [ \begin{matrix}1 & 2 \\ 0.5 & 1.0015 \end{matrix} ] , b + δb = [ \begin{matrix}3 \\ 1.4985 \end{matrix} ] .

il sistema

• (A + δA)u = b ha soluzione u = (3, 0)^{T} (perturbata matrice),

• Av = b + δb ha soluzione v = (2, 0.5)^{T} (perturbato termine noto),

• (A + δA)w = (b + δb) ha soluzione w = (5, −1)^{T} (perturbati matrice e termine noto).

Ne risulta che per possibili piccole perturbazioni alla matrice e/o al termine noto abbiamo ottenuto nei diversi casi una soluzione molto distante dall’originaria.

Con riferimento al teorema 1.10

• det(A) ≈ 0.003, det(δA) ≈ 0.0015,

• k_{2}(A) ≈ 2083.666853410356,

• abbiamo

\frac{k_{2}(A)}{1 − k_{2}(A)∥\frac{∥δA∥_{2}}{∥A∥_{2}}} ≈ 3.05387e + 04,

• \frac{∥δx∥_{2}}{∥x∥_{2}} ≈ 3.16228e + 00,

• \frac{∥δb∥_{2}}{∥b∥_{2}} ≈ 4.47214e − 04 e \frac{∥δA∥_{2}}{∥A∥_{2}} ≈ 4.47178e − 04.

La stima asserisce in effetti che

3.16228e + 00 ≤ 3.05387e + 04 · (4.47214e − 04 + 4.47178e − 04)
≈ 2.73136e + 01.  (12)

2. Metodo di eliminazione di Gauss e fattorizzazione LU

Il proposito di questa sezione è di introdurre il metodo dell’eliminazione gaussiana per risolvere sistemi lineari Ax = b dove supponiamo A invertibile, cioè tale che det(A) ≠ 0.

Per semplificare l’esposizione cominciamo con un esempio.

Esempio. Si consideri il sistema lineare

\{ \begin{matrix}3x_{1} − 2x_{2} − x_{3} = 0 \\ 6x_{1} − 2x_{2} + 2x_{3} = 6 \\ −9x_{1} + 7x_{2} + x_{3} = −1 \end{matrix}

Sappiamo che

• se moltiplichiamo ambo i membri di una stessa equazione per una costante non nulla, il nuovo problema ha le stesse soluzioni del precedente;

• se sommiamo membro a membro due equazioni distinte, il nuovo problema ha le stesse soluzioni del precedente;

• in generale se moltiplichiamo una equazione per una costante e sommiamo il risultato membro a membro a un’altra equazione, otteniamo un nuovo problema che ha le stesse soluzioni del precedente.

Diremo che due sistemi sono equivalenti, se hanno le stesse soluzioni.

• moltiplichiamo la prima riga per −2 e la sommiamo alla seconda, ottenendo

\{ \begin{matrix}3x_{1} − 2x_{2} − x_{3} = 0 \\ 0x_{1} + 2x_{2} + 4x_{3} = 6 \\ −9x_{1} + 7x_{2} + x_{3} = −1 \end{matrix}

• moltiplichiamo la prima riga per 3 e la sommiamo alla terza, ottenendo

\{ \begin{matrix}3x_{1} − 2x_{2} − x_{3} = 0 \\ 0x_{1} + 2x_{2} + 4x_{3} = 6 \\ 0x_{1} + 1x_{2} − 2x_{3} = −1 \end{matrix}

• moltiplichiamo la seconda riga per −1/2 e la sommiamo alla terza, ottenendo

\{ \begin{matrix}3x_{1} − 2x_{2} − x_{3} = 0 \\ 0x_{1} + 2x_{2} + 4x_{3} = 6 \\ 0x_{1} + 0x_{2} − 4x_{3} = −4 \end{matrix}

Il sistema finale, è facile da risolvere.

• L’ultima equazione ha una sola variabile e una sola incognita e comporta che

x_{3} = 1.

• Inseriamo questo risultato nella penultima equazione e otteniamo

2x_{2} + 4 · 1 = 6

da cui 2x_{2} = 2 e quindi x_{2} = 1.

• Inseriamo questi risultati nella prima equazione e otteniamo

3x_{1} − 2 · 1 − 1 = 0

da cui

x_{1} = 1.

Il sistema può essere scritto matricialmente come Ax = b dove

A = ( \begin{matrix}3 & −2 & −1 \\ 6 & −2 & 2 \\ −9 & 7 & 1 \end{matrix} ) , b = ( \begin{matrix}0 \\ 6 \\ 1 \end{matrix} ) ,

e abbiamo mostrato essere equivalente a

U = ( \begin{matrix}3 & −2 & −1 \\ 0 & 2 & 4 \\ 0 & 0 & −4 \end{matrix} ) , c = ( \begin{matrix}0 \\ 6 \\ −4 \end{matrix} ) .

Dal punto di vista matriciale, abbiamo generato una sequenza di sistemi lineari equivalenti, ma via via più semplici. L’ultimo è un sistema Ux = c con U triangolare superiore, ovvero tale che U_{i,j} = 0 qualora i > j, che abbiamo risolto facilmente mediante sostituzione all’indietro.

2.1. Sistemi triangolari Definizione. Una matrice A = (a_{i,j}) ∈ ℝ^{n×n} si dice

• triangolare inferiore se A_{i,j} = 0 qualora i < j;

• triangolare superiore se A_{i,j} = 0 qualora i > j.

Bibliografia

[1] K. Atkinson, Introduction to Numerical Analysis, Wiley, 1989.

[2] A. Quarteroni e F. Saleri, Elementi di calcolo numerico, Progetto Leonardo, 1999.