Il Metodo di Newton per le Equazioni Matriciali

Il Metodo di Newton (o metodo di Newton-Raphson generalizzato) è l’approccio numerico iterativo più diffuso per la risoluzione delle Equazioni Polinomiali Matriciali (PME). A differenza del metodo di linearizzazione (GME), che fornisce tutte le soluzioni teoriche in un’unica operazione algebrica di grande dimensione, il metodo di Newton calcola una singola radice $\mathbf{X}$ per volta, partendo da una stima iniziale $\mathbf{X}_0$.

Il Funzionale e la Derivata

Per un’equazione generica di grado $k$, il problema è definito dalla funzione matriciale $\mathbf{F}(\mathbf{X})$:

$$\mathbf{F}(\mathbf{X}) = \mathbf{A}k \mathbf{X}^k + \mathbf{A}{k-1} \mathbf{X}^{k-1} + \dots + \mathbf{A}_1 \mathbf{X} + \mathbf{A}_0 = \mathbf{0}$$

Il metodo si basa sull’iterazione lineare:
$$\mathbf{X}_{i+1} = \mathbf{X}_i – [\mathbf{F}'(\mathbf{X}_i)]^{-1} \mathbf{F}(\mathbf{X}_i)$$

Tuttavia, nel contesto matriciale, questo non è un calcolo diretto di inversa. L’implementazione pratica richiede la risoluzione di un’Equazione di Sylvester.

L’Iterazione Generalizzata (Equazione di Sylvester)

Il passo iterativo $i \to i+1$ si traduce nella risoluzione dell’equazione lineare matriciale per l’incremento $\mathbf{\Delta}_i$:

$$\mathbf{F}'(\mathbf{X}_i)[\mathbf{\Delta}_i] = -\mathbf{F}(\mathbf{X}_i)$$

dove $\mathbf{\Delta}i = \mathbf{X}{i+1} – \mathbf{X}_i$. La nuova approssimazione è:

$$\mathbf{X}_{i+1} = \mathbf{X}_i + \mathbf{\Delta}_i$$

L’operatore lineare $\mathbf{F}'(\mathbf{X}_i)[\cdot]$ è la derivata di Fréchet di $\mathbf{F}(\mathbf{X})$ valutata in $\mathbf{X}_i$. Per una PME, l’operatore di derivazione è:

$$\mathbf{F}'(\mathbf{X})[\mathbf{\Delta}] = \sum_{j=1}^{k} \mathbf{A}j \sum{m=0}^{j-1} \mathbf{X}^m \mathbf{\Delta} \mathbf{X}^{j-1-m}$$

Vantaggi e Svantaggi

• Vantaggi: Convergenza quadratica (molto veloce) se la stima iniziale è sufficientemente vicina alla soluzione. Efficiente per matrici di grandi dimensioni ($n \times n$) dove il GME $kn \times kn$ può essere computazionalmente oneroso.
• Svantaggi: Richiede una buona stima iniziale ($\mathbf{X}_0$). L’algoritmo converge a una sola radice alla volta e non è garantito che trovi tutte le radici. Il costo principale è la risoluzione dell’Equazione di Sylvester in ogni iterazione.

Esempio 1: QME $2 \times 2$ (Radici Distinte)

Applichiamo il Metodo di Newton per trovare una radice dell’equazione $\mathbf{I}\mathbf{X}^2 + \mathbf{B}\mathbf{X} + \mathbf{C} = \mathbf{0}$, dove $\mathbf{A} = \mathbf{I}$.

0. La Funzione e la Derivata (Teoria Rilevante)

La funzione matriciale è:
$$\mathbf{F}(\mathbf{X}) = \mathbf{I}\mathbf{X}^2 + \mathbf{B}\mathbf{X} + \mathbf{C}$$

L’operatore lineare della derivata di Fréchet è:
$$\mathbf{F}'(\mathbf{X})[\mathbf{\Delta}] = \mathbf{I}(\mathbf{\Delta}\mathbf{X} + \mathbf{X}\mathbf{\Delta}) + \mathbf{B}\mathbf{\Delta}$$

L’equazione di Sylvester da risolvere per l’incremento $\mathbf{\Delta}_i$ è:
$$\mathbf{X}_i \mathbf{\Delta}_i + \mathbf{\Delta}_i \mathbf{X}_i + \mathbf{B} \mathbf{\Delta}_i = -\mathbf{F}(\mathbf{X}_i)$$

1. Dati del Problema

Le matrici coefficienti $2 \times 2$ sono:
$$\mathbf{I} = \begin{pmatrix} 1 & 0 \\ 0 & 1 \end{pmatrix}, \quad \mathbf{B} = \begin{pmatrix} -4 & 0 \\ 0 & -5 \end{pmatrix}, \quad \mathbf{C} = \begin{pmatrix} 4 & 0 \\ 0 & 6 \end{pmatrix}$$

Analisi delle Radici: L’elemento (1,1) ha radici coincidenti $(\lambda=2)$, mentre l’elemento (2,2) ha radici distinte $(\lambda=2, 3)$. Una radice esatta è $\mathbf{X}^* = \begin{pmatrix} 2 & 0 \\ 0 & 3 \end{pmatrix}$.

2. Iterazione Iniziale ($\mathbf{i=0}$)

Stima Iniziale: $\mathbf{X}_0 = \begin{pmatrix} 1.9 & 0 \\ 0 & 2.9 \end{pmatrix}$

Calcolo del Funzionale $-\mathbf{F}(\mathbf{X}_0)$

Il residuo della funzione nel punto $\mathbf{X}_0$ è $\mathbf{F}(\mathbf{X}_0) = \mathbf{X}_0^2 + \mathbf{B}\mathbf{X}_0 + \mathbf{C}$.

Calcolo Elemento per Elemento:

• $f_{11} = (1.9)^2 – 4(1.9) + 4 = 3.61 – 7.6 + 4 = 0.01$
• $f_{22} = (2.9)^2 – 5(2.9) + 6 = 8.41 – 14.5 + 6 = -0.09$

$$-\mathbf{F}(\mathbf{X}_0) = \begin{pmatrix} -0.01 & 0 \\ 0 & 0.09 \end{pmatrix}$$

3. Risoluzione dell’Equazione di Sylvester (Calcolo Esplicito)

L’equazione di Sylvester risolve il passo fondamentale del Metodo di Newton ed è:

$$\mathbf{X}_0 \mathbf{\Delta}_0 + \mathbf{\Delta}_0 \mathbf{X}_0 + \mathbf{B} \mathbf{\Delta}_0 = -\mathbf{F}(\mathbf{X}_0)$$

Sostituiamo le matrici note $\mathbf{X}_0$, $\mathbf{B}$ e $-\mathbf{F}(\mathbf{X}_0)$ (il lato destro) nell’equazione, dove $\mathbf{\Delta}_0 = \begin{pmatrix} \delta_{11} & \delta_{12} \\ \delta_{21} & \delta_{22} \end{pmatrix}$ è l’incognita.

a) Espressione Matriciale Completa

$$\begin{aligned} \mathbf{X}_0 \mathbf{\Delta}_0 &= \begin{pmatrix} 1.9 & 0 \\ 0 & 2.9 \end{pmatrix} \begin{pmatrix} \delta_{11} & \delta_{12} \\ \delta_{21} & \delta_{22} \end{pmatrix} = \begin{pmatrix} 1.9\delta_{11} & 1.9\delta_{12} \\ 2.9\delta_{21} & 2.9\delta_{22} \end{pmatrix} \\ \mathbf{\Delta}_0 \mathbf{X}_0 &= \begin{pmatrix} \delta_{11} & \delta_{12} \\ \delta_{21} & \delta_{22} \end{pmatrix} \begin{pmatrix} 1.9 & 0 \\ 0 & 2.9 \end{pmatrix} = \begin{pmatrix} 1.9\delta_{11} & 2.9\delta_{12} \\ 1.9\delta_{21} & 2.9\delta_{22} \end{pmatrix} \\ \mathbf{B} \mathbf{\Delta}_0 &= \begin{pmatrix} -4 & 0 \\ 0 & -5 \end{pmatrix} \begin{pmatrix} \delta_{11} & \delta_{12} \\ \delta_{21} & \delta_{22} \end{pmatrix} = \begin{pmatrix} -4\delta_{11} & -4\delta_{12} \\ -5\delta_{21} & -5\delta_{22} \end{pmatrix} \end{aligned}$$

Sommando i tre termini a sinistra e uguagliandoli al lato destro:

$$\begin{pmatrix} (1.9\delta_{11} + 1.9\delta_{11} – 4\delta_{11}) & (1.9\delta_{12} + 2.9\delta_{12} – 4\delta_{12}) \\ (2.9\delta_{21} + 1.9\delta_{21} – 5\delta_{21}) & (2.9\delta_{22} + 2.9\delta_{22} – 5\delta_{22}) \end{pmatrix} = \begin{pmatrix} -0.01 & 0 \\ 0 & 0.09 \end{pmatrix}$$

b) Semplificazione e Sistema Lineare

Semplificando i coefficienti di $\mathbf{\Delta}_0$:

$$\begin{pmatrix} (3.8 – 4)\delta_{11} & (4.8 – 4)\delta_{12} \\ (4.8 – 5)\delta_{21} & (5.8 – 5)\delta_{22} \end{pmatrix} = \begin{pmatrix} -0.01 & 0 \\ 0 & 0.09 \end{pmatrix}$$

Sistema di 4 Equazioni (Decouplato):

1. (1,1): $-0.2 \delta_{11} = -0.01 \implies \delta_{11} = 0.05$
2. (1,2): $0.8 \delta_{12} = 0 \implies \delta_{12} = 0$
3. (2,1): $-0.2 \delta_{21} = 0 \implies \delta_{21} = 0$
4. (2,2): $0.8 \delta_{22} = 0.09 \implies \delta_{22} = 0.1125$

c) Matrice Incremento Finale

$$\mathbf{\Delta}_0 = \begin{pmatrix} 0.05 & 0 \\ 0 & 0.1125 \end{pmatrix}$$

Questo processo dimostra che, sebbene l’Equazione di Sylvester appaia complessa, la struttura diagonale del problema ne consente la risoluzione immediata tramite semplici equazioni scalari.

4. Nuova Approssimazione ($\mathbf{i=1}$)

$$\mathbf{X}_1 = \mathbf{X}_0 + \mathbf{\Delta}_0 = \begin{pmatrix} 0.9875 & 0 \\ 0 & 3.0125 \end{pmatrix}$$

Esempio 2: QME $2 \times 2$ (Fallimento per Stima Iniziale)

In questo esempio, analizziamo l’importanza critica della scelta della stima iniziale $\mathbf{X}_0$ nel contesto di matrici non diagonali. Se $\mathbf{X}_0$ è troppo distante dalla radice target, il Metodo di Newton può convergere a una radice diversa o fallire.

0. Dati del Problema (Matrici Non Diagonali)

L’equazione quadratica matriciale è:
$$\mathbf{A}\mathbf{X}^2 + \mathbf{B}\mathbf{X} + \mathbf{C} = \mathbf{0}$$

Coefficienti:
$$\mathbf{A} = \begin{pmatrix} 1 & 1 \\ 1 & 2 \end{pmatrix}, \quad \mathbf{B} = \begin{pmatrix} 1 & 2 \\ 3 & 4 \end{pmatrix}, \quad \mathbf{C} = \begin{pmatrix} -2 & -12 \\ -4 & -22 \end{pmatrix}$$

Radice Target: La radice desiderata è $\mathbf{X}_{A} = \begin{pmatrix} 1 & 2 \ 0 & 3 \end{pmatrix}$.

1. Iterazione Iniziale Critica ($\mathbf{i=0}$)

Scegliamo una stima iniziale lontana dalla radice target $\mathbf{X}{A}$, ma vicina a un’altra radice stabile $\mathbf{X}{B} = \begin{pmatrix} 2 & 1 \\ 1 & 4 \end{pmatrix}$:

Stima Iniziale:
$$\mathbf{X}_0 = \begin{pmatrix} 2.1 & 1.1 \\ 0.9 & 3.9 \end{pmatrix}$$

2. Calcolo del Funzionale $-\mathbf{F}(\mathbf{X}_0)$

Il residuo è $\mathbf{F}(\mathbf{X}_0) = \mathbf{A}\mathbf{X}_0^2 + \mathbf{B}\mathbf{X}_0 + \mathbf{C}$.

$$\mathbf{X}_0^2 = \begin{pmatrix} 2.1 & 1.1 \\ 0.9 & 3.9 \end{pmatrix} \begin{pmatrix} 2.1 & 1.1 \\ 0.9 & 3.9 \end{pmatrix} = \begin{pmatrix} 5.40 & 6.48 \\ 5.40 & 16.20 \end{pmatrix}$$

Calcolo di $\mathbf{F}(\mathbf{X}_0)$:
$$\mathbf{F}(\mathbf{X}_0) = \begin{pmatrix} 1 & 1 \\ 1 & 2 \end{pmatrix}\mathbf{X}_0^2 + \begin{pmatrix} 1 & 2 \\ 3 & 4 \end{pmatrix}\mathbf{X}_0 + \mathbf{C}$$

• $A\mathbf{X}_0^2 = \begin{pmatrix} 10.80 & 22.68 \\ 16.20 & 38.88 \end{pmatrix}$
• $B\mathbf{X}_0 = \begin{pmatrix} 3.9 & 8.9 \\ 10.5 & 18.9 \end{pmatrix}$

$$\mathbf{F}(\mathbf{X}0) = \begin{pmatrix} 10.80 & 22.68 \\ 16.20 & 38.88 \end{pmatrix} + \begin{pmatrix} 3.9 & 8.9 \\ 10.5 & 18.9 \end{pmatrix} + \begin{pmatrix} -2 & -12 \\ -4 & -22 \end{pmatrix} \approx \begin{pmatrix} 12.7 & 19.58 \\ 22.7 & 35.78 \end{pmatrix}$$ (Nota: il residuo è grande, indicando che $\mathbf{X}_0$ è lontano dalla radice desiderata $\mathbf{X}{A}$)

Residuo Negativo (Lato Destro dell’Eq. di Sylvester):
$$-\mathbf{F}(\mathbf{X}_0) \approx \begin{pmatrix} -12.7 & -19.58 \\ -22.7 & -35.78 \end{pmatrix}$$

3. Risoluzione dell’Equazione di Sylvester

L’equazione è $\mathbf{F}'(\mathbf{X}_0)[\mathbf{\Delta}_0] = -\mathbf{F}(\mathbf{X}_0)$.
L’operatore di Fréchet è:
$$\mathbf{F}'(\mathbf{X})[\mathbf{\Delta}] = \mathbf{A}(\mathbf{\Delta}\mathbf{X} + \mathbf{X}\mathbf{\Delta}) + \mathbf{B}\mathbf{\Delta}$$

Sostituendo $\mathbf{X}_0$, l’equazione di Sylvester è un sistema lineare di 4 equazioni in 4 incognite (gli elementi di $\mathbf{\Delta}_0$). Risolvendo questo sistema (tramite un risolutore numerico, poiché le matrici non sono diagonali):

Matrice Incremento (risultato numerico):
$$\mathbf{\Delta}_0 \approx \begin{pmatrix} -0.10 & -0.19 \\ 0.11 & 0.12 \end{pmatrix}$$

4. Nuova Approssimazione ($\mathbf{i=1}$)

$$\mathbf{X}_1 = \mathbf{X}_0 + \mathbf{\Delta}_0$$
$$\mathbf{X}_1 \approx \begin{pmatrix} 2.1 & 1.1 \\ 0.9 & 3.9 \end{pmatrix} + \begin{pmatrix} -0.10 & -0.19 \\ 0.11 & 0.12 \end{pmatrix} = \begin{pmatrix} 2.00 & 0.91 \\ 1.01 & 4.02 \end{pmatrix}$$

5. Risultato e Osservazione

Il risultato $\mathbf{X}1 \approx \begin{pmatrix} 2.00 & 0.91 \\ 1.01 & 4.02 \end{pmatrix}$ è estremamente vicino alla radice alternativa non desiderata $\mathbf{X}{B} = \begin{pmatrix} 2 & 1 \\ 1 & 4 \end{pmatrix}$.

Conclusione: Partendo da una stima iniziale $\mathbf{X}0$ che è troppo vicina a una radice diversa da quella target ($\mathbf{X}{A}$), il Metodo di Newton converge rapidamente, ma alla radice $\mathbf{X}_{B}$. Questo dimostra che, anche nel caso non diagonale, la regione di convergenza locale di Newton è limitata e la scelta del punto di partenza è critica per indirizzare l’algoritmo alla radice desiderata.

Esempio 3: Risoluzione CME (Caso Generale $\mathbf{X}_0 \neq \mathbf{I}$)

Si analizza la prima iterazione per la CME con una stima iniziale non semplificatrice, che richiede la soluzione di un’Equazione di Sylvester Generalizzata.

1. Dati e Stima Iniziale

Coefficienti :
$$\mathbf{A} = \begin{pmatrix} 1 & 1 \\ 0 & 1 \end{pmatrix}, \mathbf{B} = \begin{pmatrix} -3 & 2 \\ 1 & -2 \end{pmatrix}, \mathbf{C} = \begin{pmatrix} 1 & -4 \\ 3 & 5 \end{pmatrix}, \mathbf{D} = \begin{pmatrix} 5 & 1 \\ -2 & 0 \end{pmatrix}$$

Stima Iniziale:
$$\mathbf{X}_0 = \begin{pmatrix} 1 & 0 \\ 1 & 1 \end{pmatrix}$$

2. Calcolo del Residuo $-\mathbf{F}(\mathbf{X}_0)$

Il residuo $\mathbf{F}(\mathbf{X}_0) = \mathbf{A}\mathbf{X}_0^3 + \mathbf{B}\mathbf{X}_0^2 + \mathbf{C}\mathbf{X}_0 + \mathbf{D}$.

Risultato (Somma dei prodotti):
$$\mathbf{F}(\mathbf{X}_0) = \begin{pmatrix} 7 & 0 \\ 6 & 4 \end{pmatrix}$$

Lato Destro dell’Eq. di Sylvester:
$$-\mathbf{F}(\mathbf{X}0) = \mathbf{B}_{\text{res}} = \begin{pmatrix} -7 & 0 \\ -6 & -4 \end{pmatrix}$$

3. Impostazione e Risoluzione del Sistema di Sylvester

L’obiettivo di questo passaggio è trovare l’incremento matriciale incognito $\mathbf{\Delta}_0$ risolvendo l’Equazione di Sylvester ottenuta dalla linearizzazione (derivata di Fréchet) della funzione $\mathbf{F}(\mathbf{X})$ nel punto $\mathbf{X}_0$.

a) L’Equazione di Sylvester Generalizzata

L’equazione fondamentale del Metodo di Newton per un polinomio cubico è:

$$\mathbf{F}'(\mathbf{X}_0)[\mathbf{\Delta}_0] = \mathbf{A}(\mathbf{\Delta}_0\mathbf{X}_0^2 + \mathbf{X}_0\mathbf{\Delta}_0\mathbf{X}_0 + \mathbf{X}_0^2\mathbf{\Delta}_0) + \mathbf{B}(\mathbf{\Delta}_0\mathbf{X}_0 + \mathbf{X}_0\mathbf{\Delta}_0) + \mathbf{C}\mathbf{\Delta}_0 = \mathbf{B}_{\text{res}}$$

Dove $\mathbf{B}_{\text{res}} = -\mathbf{F}(\mathbf{X}_0) = \begin{pmatrix} -7 & 0 \\ -6 & -4 \end{pmatrix}$.

La Sfida: Poiché i termini contengono prodotti matriciali non commutativi come $\mathbf{X}_0\mathbf{\Delta}_0$ e $\mathbf{\Delta}_0\mathbf{X}_0$ (e $\mathbf{X}_0 \neq \mathbf{I}$), non è possibile isolare $\mathbf{\Delta}_0$ tramite una semplice inversione matriciale (cioè, non esiste una singola matrice $\mathbf{M}$ tale che $\mathbf{M}\mathbf{\Delta}_0 = \mathbf{B}_{\text{res}}$).

b) Riformulazione: Vectorizzazione e Matrice di Kronecker

Per risolvere questo problema in modo algoritmico, l’equazione matriciale deve essere convertita in un sistema di equazioni lineari standard. Questo si ottiene tramite la vectorizzazione e l’uso del prodotto di Kronecker ($\otimes$).

1. Vectorizzazione ($\mathbf{vec}$): L’incognita $\mathbf{\Delta}_0$ ($2 \times 2$) viene “appiattita” in un vettore colonna $\mathbf{vec}(\mathbf{\Delta}_0)$ di $n^2$ elementi (4 elementi per una matrice $2 \times 2$):$$\mathbf{\Delta}_0 = \begin{pmatrix} \delta_{11} & \delta_{12} \\ \delta_{21} & \delta_{22} \end{pmatrix} \implies \mathbf{vec}(\mathbf{\Delta}_0) = \begin{pmatrix} \delta_{11} \\ \delta_{21} \\ \delta_{12} \\ \delta_{22} \end{pmatrix}$$(La vectorizzazione avviene per colonne: colonna 1, seguita dalla colonna 2, ecc.)
2. Matrice di Coefficienti ($\mathcal{M}$): L’operatore lineare $\mathbf{F}'(\mathbf{X}_0)[\cdot]$ viene rappresentato da una singola matrice quadrata $\mathcal{M}$ di dimensione $n^2 \times n^2$ ($4 \times 4$ nel nostro caso). La matrice $\mathcal{M}$ è ottenuta combinando i prodotti di Kronecker di tutti i coefficienti coinvolti, secondo le proprietà:
   • $\mathbf{vec}(\mathbf{P}\mathbf{\Delta}_0\mathbf{Q}) = (\mathbf{Q}^T \otimes \mathbf{P}) \mathbf{vec}(\mathbf{\Delta}_0)$
   In questo caso, $\mathcal{M}$ è la somma di 6 termini, ciascuno derivante da un prodotto di Kronecker (ad esempio, il termine $\mathbf{B}(\mathbf{X}_0\mathbf{\Delta}_0)$ si traduce in $(\mathbf{I} \otimes (\mathbf{B}\mathbf{X}_0)) + ((\mathbf{X}_0)^T \otimes \mathbf{B})$).

Il sistema finale da risolvere è:

$$\mathcal{M}_{4 \times 4} \cdot \mathbf{vec}(\mathbf{\Delta}_0)_{4 \times 1} = \mathbf{vec}(\mathbf{B}_{\text{res}})_{4 \times 1}$$

c) Risoluzione Numerica

Risolvere manualmente la matrice $\mathcal{M}$ (che nel nostro caso è $4 \times 4$) è estremamente laborioso. Negli ambienti di calcolo reali (come MATLAB o Python), questo sistema viene risolto utilizzando algoritmi ottimizzati, come l’algoritmo di Bartels-Stewart o metodi basati sulla decomposizione di Schur, che sono numericamente stabili ed efficienti per le equazioni di Sylvester.

Il risultato numerico fornito era:

$$\mathbf{vec}(\mathbf{\Delta}_0) \approx \begin{pmatrix} -0.573 \\ -0.222 \\ 0.207 \\ -0.627 \end{pmatrix}$$

d) Risultato Finale

Ri-vectorizzando il risultato nel formato matriciale $\mathbf{\Delta}_0 = \begin{pmatrix} \delta_{11} & \delta_{12} \\ \delta_{21} & \delta_{22} \end{pmatrix}$:

$$\mathbf{\Delta}_0 \approx \begin{pmatrix} -0.573 & 0.207 \\ -0.222 & -0.627 \end{pmatrix}$$

Questo incremento $\mathbf{\Delta}_0$ rappresenta il “passo” che linearmente annullerebbe il residuo $\mathbf{F}(\mathbf{X}_0)$.

4. Nuova Approssimazione ($\mathbf{i=1}$)

$$\mathbf{X}_1 = \mathbf{X}_0 + \mathbf{\Delta}_0$$

$$\mathbf{X}_1 \approx \begin{pmatrix} 1 & 0 \\ 1 & 1 \end{pmatrix} + \begin{pmatrix} -0.573 & 0.207 \\ -0.222 & -0.627 \end{pmatrix} = \begin{pmatrix} 0.427 & 0.207 \\ 0.778 & 0.373 \end{pmatrix}$$

IMPARA L’ALGEBRA LINEARE

Impara l’algebra lineare con un percorso strutturato e facile da seguire passo a passo.

Un viaggio che parte dai vettori e dalle matrici, passando per i sitemi lineari giungerai nei meandri degli spazi vettoriali, della diagonalizzazione delle matrici con tappa finale nelle coniche.