Algoritmo di Gauss e fattorizzazione LU
Supponiamo di voler fattorizzare una matrice $ntimesn$ in forma LU (A=LU, L speciale triangolare inferiore, U triangolare superiore). Sappiamo che:
*) non sempre è possibile, ma che sicuramente esiste una matrice di permutazione $Pi$ tale che $PiA$ si fattorizza.
L'algoritmo di eliminazione di Gauss applicato ad una matrice quadrata produce (quando applicabile) una successione di mat. elementari d Gauss tali che $E^((n-1))...E^((1))A=U$, e il prodotto delle $E^((i))$ è una matrice speciale triangolare inferiore. Quindi l'algoritmo di Gauss, quando applicabile, è (anche) un metodo per calcolare la fattorizzazione LU.
L'algoritmo di Gauss diventa sicuramente applicabile se ammettiamo la possibilità di scambiare delle righe ad ogni passo. Quindi il risultato sarà una cosa tipo $E^((n-1))Pi^((n-1))E^((n-2))...E^((1))Pi^((1))A=U$.
Domanda: c'entra qualcosa con il risultato segnato con *)? Penso proprio di sì. Infatti vorrei dimostrare che la produttoria di matrici $E^((n-1))Pi^((n-1))E^((n-2))...E^((1))Pi^((1))$ generata dall'algoritmo di Gauss con scambi sia uguale ad una cosa come $M^((n-1))M^((n-2))...M^((1))*Sigma^((n-1))...Sigma^((1))$ dove le $M$ sono matrici di Gauss e le $Sigma$ sono matrici di permutazione. (Questo significherebbe che, in linea teorica, sarebbe possibile effettuare gli scambi di riga all'inizio dell'algoritmo in un'unica soluzione, e poi applicare le trasformazioni di Gauss senza scambi).
Potrebbe funzionare o sto perdendo tempo?
*) non sempre è possibile, ma che sicuramente esiste una matrice di permutazione $Pi$ tale che $PiA$ si fattorizza.
L'algoritmo di eliminazione di Gauss applicato ad una matrice quadrata produce (quando applicabile) una successione di mat. elementari d Gauss tali che $E^((n-1))...E^((1))A=U$, e il prodotto delle $E^((i))$ è una matrice speciale triangolare inferiore. Quindi l'algoritmo di Gauss, quando applicabile, è (anche) un metodo per calcolare la fattorizzazione LU.
L'algoritmo di Gauss diventa sicuramente applicabile se ammettiamo la possibilità di scambiare delle righe ad ogni passo. Quindi il risultato sarà una cosa tipo $E^((n-1))Pi^((n-1))E^((n-2))...E^((1))Pi^((1))A=U$.
Domanda: c'entra qualcosa con il risultato segnato con *)? Penso proprio di sì. Infatti vorrei dimostrare che la produttoria di matrici $E^((n-1))Pi^((n-1))E^((n-2))...E^((1))Pi^((1))$ generata dall'algoritmo di Gauss con scambi sia uguale ad una cosa come $M^((n-1))M^((n-2))...M^((1))*Sigma^((n-1))...Sigma^((1))$ dove le $M$ sono matrici di Gauss e le $Sigma$ sono matrici di permutazione. (Questo significherebbe che, in linea teorica, sarebbe possibile effettuare gli scambi di riga all'inizio dell'algoritmo in un'unica soluzione, e poi applicare le trasformazioni di Gauss senza scambi).
Potrebbe funzionare o sto perdendo tempo?
Risposte
$E^((n-1))Pi^((n-1))E^((n-2))...E^((1))Pi^((1))=>M^((n-1))M^((n-2))...M^((1))*Sigma^((n-1))...Sigma^((1))$ lo puoi ottenere, basta che ricordi che le matrici di permutazione scambiano le righe e quindi applicate due volte torni allo scambio iniziale, cioè $(Pi^(j))^2=I_d$. Basta giocare con questo fatto.
infatti per esempio il caso con due matrici avresti $M_2P_2M_1P_1A=M_2P_2M_1P_2P_2P_1A=M_2hatM_1P_2P_1A=U
con $hatM_1=P_2M_1P_2$ che non è una matrice come quella di prima, però è dello stesso tipo e quindi va bene
, inoltre $P_2P_1=P$ è la matrice di permutazione tale che $PA=LU$
in questo modo se hai n matrici ottiene $hatM_nP_n...P_1A=U$ con $hatM_n=M_nhatM_(n-1)hathatM_(n-2)...hathathathatM_1$ con
$hatM_(n_1)=P_nM_(n-1)P_n
$hathatM_(n-2)=P_nP_(n-1)M_(n-2)P_nP_(n-1)
...
$hathathathatM_1=P_n...P_2M_1P_n...P_2
quindi ottieni infine che $L=hatM_n^(-1)$ quindi volendo puoi fare quello che dici, però sarebbe una follia dal punto di vista computazionale se ci pensi...
nel senso in questo modo al posto di scambiare le righe ad ogni passo devi fare due volte una moltiplicazione tra matrici... e devi anche pensare che le matrici $M_j$ si formano man mano che va avanti la composizione, quindi la P si forma man mano...
in conclusione magari si può, ma
a) non avrebbe un vantaggio quando costruisci l'algoritmo
b) devi riuscire a capire un modo che ti permetta di sapere in anticipo come si scambino le righe mentre si formano le matrici di Gauss.
anche se normalemnte la differenza non è abissale in quanto costruisci separatamente la P e la L in un programma generico.
spero di essere stato utile e chiaro.
infatti per esempio il caso con due matrici avresti $M_2P_2M_1P_1A=M_2P_2M_1P_2P_2P_1A=M_2hatM_1P_2P_1A=U
con $hatM_1=P_2M_1P_2$ che non è una matrice come quella di prima, però è dello stesso tipo e quindi va bene
, inoltre $P_2P_1=P$ è la matrice di permutazione tale che $PA=LU$in questo modo se hai n matrici ottiene $hatM_nP_n...P_1A=U$ con $hatM_n=M_nhatM_(n-1)hathatM_(n-2)...hathathathatM_1$ con
$hatM_(n_1)=P_nM_(n-1)P_n
$hathatM_(n-2)=P_nP_(n-1)M_(n-2)P_nP_(n-1)
...
$hathathathatM_1=P_n...P_2M_1P_n...P_2
quindi ottieni infine che $L=hatM_n^(-1)$ quindi volendo puoi fare quello che dici, però sarebbe una follia dal punto di vista computazionale se ci pensi...
nel senso in questo modo al posto di scambiare le righe ad ogni passo devi fare due volte una moltiplicazione tra matrici... e devi anche pensare che le matrici $M_j$ si formano man mano che va avanti la composizione, quindi la P si forma man mano...
in conclusione magari si può, ma
a) non avrebbe un vantaggio quando costruisci l'algoritmo
b) devi riuscire a capire un modo che ti permetta di sapere in anticipo come si scambino le righe mentre si formano le matrici di Gauss.
anche se normalemnte la differenza non è abissale in quanto costruisci separatamente la P e la L in un programma generico.
spero di essere stato utile e chiaro.
Certo. Hai ragione: dal punto di vista computazionale è una fesseria.
Però forse riusciamo a ricavarci una dimostrazione più facile per il teorema di esistenza della fattorizzazione LU a meno di matrici di permutazione (quello che avevo segnato con *). Infatti:
sul fatto che esiste una successione ${E^((i))\ |\ i=1..n}$ di matrici di Gauss e una ${Pi^((j))|\ j=1..n}$ di matrici di permutazione (che di fatto sono solo scambi singoli, non permutazioni più complicate) tale che $E^((n-1))Pi^((n-1))E^((n-2))...E^((1))Pi^((1))A=U$ triang. superiore non ci piove.
se dimostriamo che possiamo "separare" le matrici di Gauss da quelle di permutazione -come hai fatto nel precedente post- abbiamo finito: abbiamo trovato una mat. di permutazione $Pi$ tale che $PiA$ si fattorizza in LU. (E' una questione solo teorica, per intenderci.)
Resta da dimostrare per bene quest'ultimo punto. Tu l'hai fatto ma non mi è chiarissimo il meccanismo, per un solo motivo: come fai a dire che $hatM_1$ è ancora una matrice di Gauss?
Però forse riusciamo a ricavarci una dimostrazione più facile per il teorema di esistenza della fattorizzazione LU a meno di matrici di permutazione (quello che avevo segnato con *). Infatti:
sul fatto che esiste una successione ${E^((i))\ |\ i=1..n}$ di matrici di Gauss e una ${Pi^((j))|\ j=1..n}$ di matrici di permutazione (che di fatto sono solo scambi singoli, non permutazioni più complicate) tale che $E^((n-1))Pi^((n-1))E^((n-2))...E^((1))Pi^((1))A=U$ triang. superiore non ci piove.
se dimostriamo che possiamo "separare" le matrici di Gauss da quelle di permutazione -come hai fatto nel precedente post- abbiamo finito: abbiamo trovato una mat. di permutazione $Pi$ tale che $PiA$ si fattorizza in LU. (E' una questione solo teorica, per intenderci.)
Resta da dimostrare per bene quest'ultimo punto. Tu l'hai fatto ma non mi è chiarissimo il meccanismo, per un solo motivo: come fai a dire che $hatM_1$ è ancora una matrice di Gauss?
perchè viene moltiplicata per matrici che scambiano solo le righe e quindi M è composta solo dai moltiplicatori come nella matrice di Gauss.
E allora non va bene secondo me. Ad esempio, la matrice $[[1,0,0],[0,1,0],[0,1,1]]$ è di Gauss. Prendiamo la matrice $Pi=[[0,0,1],[0,1,0],[1,0,0]]$, che è uno scambio (moltiplicata a sinistra, scambia la prima e la terza riga). Infatti $Pi*Pi=I$.
Però prova a fare $Pi*[[1,0,0],[0,1,0],[0,1,1]]*Pi$: ottieni $[[1,1,0],[0,1,0],[0,0,1]]$ che non è una matrice di Gauss. Non è nemmeno triangolare inferiore, quindi manda gambe all'aria pure la fattorizzazione LU. Forse la proposizione è falsa, sai?
Però prova a fare $Pi*[[1,0,0],[0,1,0],[0,1,1]]*Pi$: ottieni $[[1,1,0],[0,1,0],[0,0,1]]$ che non è una matrice di Gauss. Non è nemmeno triangolare inferiore, quindi manda gambe all'aria pure la fattorizzazione LU. Forse la proposizione è falsa, sai?
effettivamente hai ragione, facendo PMP ci si trova ad avere una triangolare superiore... mmmm mi era sfuggito... devo riflettere per un metodo alternativo. Separare dovresti riuscire a separare in quanto dalla sequenza iniziale di P ed M mescolate devi isolare le M e le P, così facendo riesci poi a ricavare L...
ci penso su
ci penso su
ci ho pensato... allora abbiamo che una generica matrice $M_j*P_j$ è una matrice di Gauss con al più le righe scambiate. Queste creano la successione di matrici richiesta.
condieriamo $M_n*P_n*M_(n-1)*P_(n-1)...M_2*P_2*M_1*P_1*A=U$ ricordandoci come prima che $P_j^2=I_d$ possiamo scrivere:
$M_n*P_n*M_(n-1)*P_(n-1)...M_2*P_2*M_1*P_2*P_3...P_(n-1)*P_n*P_n*P_(n-1)...P_3*P_2*P_1*A=U
poniamo $P_2P_3...P_(n_1)P_n=hatP$ . Abbiamo ottenuto la seguente ugualianza:
$M_n*P_n*M_(n-1)*P_(n-1)...M_2*P_2*M_1*hatP*P_n*P_(n-1)...P_3*P_2*P_1*A=U
ricordando quello detto all'inizio, cioè che $M_j*P_j=hatM_j$ è una matrice di Gauss con al più le righe scambiate, lo è in particolare anche $M_1*hatP$. Posto $P_n*P_(n-1)...P_3*P_2*P_1=Pi$ otteniamo
$hatM_m*hatM_(n-1)...hatM_1*Pi*A=U$ come speravamo.
Più convincente ora?
condieriamo $M_n*P_n*M_(n-1)*P_(n-1)...M_2*P_2*M_1*P_1*A=U$ ricordandoci come prima che $P_j^2=I_d$ possiamo scrivere:
$M_n*P_n*M_(n-1)*P_(n-1)...M_2*P_2*M_1*P_2*P_3...P_(n-1)*P_n*P_n*P_(n-1)...P_3*P_2*P_1*A=U
poniamo $P_2P_3...P_(n_1)P_n=hatP$ . Abbiamo ottenuto la seguente ugualianza:
$M_n*P_n*M_(n-1)*P_(n-1)...M_2*P_2*M_1*hatP*P_n*P_(n-1)...P_3*P_2*P_1*A=U
ricordando quello detto all'inizio, cioè che $M_j*P_j=hatM_j$ è una matrice di Gauss con al più le righe scambiate, lo è in particolare anche $M_1*hatP$. Posto $P_n*P_(n-1)...P_3*P_2*P_1=Pi$ otteniamo
$hatM_m*hatM_(n-1)...hatM_1*Pi*A=U$ come speravamo.
Più convincente ora?
ehi, grazie fu! io veramente ci avevo messo una croce sopra, ma questo risultato adesso mi serve perché vorrei mettere a punto un algoritmo di fattorizzazione LU con variante del pivot (parziale).
In pratica dovrei determinare esplicitamente le $L, U$ facendo però scambi di riga in maniera tale da contenerne la crescita degli elementi. Il problema è che così non ottengo la fattorizzazione LU della matrice originaria, ma quella di qualcos'altro che vorrei riuscire a determinare. Precisamente, usando la tua scrittura, ottengo $\hat{M_n}...\hat{M_1}PiA=U$ e quindi la fattorizzazione $PiA=(\hat{M_n}...\hat{M_1})^(-1)U$. Il fatto che ci sia $Pi$ non è un problema perché io so come è fatta (sarà il prodotto di tutti gli scambi fatti durante il procedimento).
Invece il problema è che non sono convinto del fatto che le $\hat{M_i}$ sono triangolari inferiori. Perché se tu a una matrice triangolare inferiore permuti le righe ottieni qualcos'altro... Devo pensarci un po', appena ho tempo.
In pratica dovrei determinare esplicitamente le $L, U$ facendo però scambi di riga in maniera tale da contenerne la crescita degli elementi. Il problema è che così non ottengo la fattorizzazione LU della matrice originaria, ma quella di qualcos'altro che vorrei riuscire a determinare. Precisamente, usando la tua scrittura, ottengo $\hat{M_n}...\hat{M_1}PiA=U$ e quindi la fattorizzazione $PiA=(\hat{M_n}...\hat{M_1})^(-1)U$. Il fatto che ci sia $Pi$ non è un problema perché io so come è fatta (sarà il prodotto di tutti gli scambi fatti durante il procedimento).
Invece il problema è che non sono convinto del fatto che le $\hat{M_i}$ sono triangolari inferiori. Perché se tu a una matrice triangolare inferiore permuti le righe ottieni qualcos'altro... Devo pensarci un po', appena ho tempo.
@fu^2: Ho trovato una dispensa online, all'indirizzo http://docenti.ing.unipi.it/~d3253/libro/sistemi.pdf .
A pagina 6 del pdf (pag. 54 totale) leggo:
Quindi la proposizione è vera. Solo che dal tuo post non riesco a capire, oltre al fatto delle matrici triangolari inferiori che dicevo prima, quale sia la matrice di permutazione P: e in fondo, questo è proprio quello che mi serve.
A pagina 6 del pdf (pag. 54 totale) leggo:
Nel caso di una matrice A qualunque si può dimostrare che l’algoritmo
di Gauss con l’eventuale uso del pivoting parziale conduce ancora ad una
fattorizzazione della forma
PA = LpRp (3.15)
dove P è una matrice di permutazione definita dagli scambi di righe richiesti
dall’algoritmo, Rp è triangolare superiore ed Lp è triangolare inferiore con
elementi diagonali unitari.
Quindi la proposizione è vera. Solo che dal tuo post non riesco a capire, oltre al fatto delle matrici triangolari inferiori che dicevo prima, quale sia la matrice di permutazione P: e in fondo, questo è proprio quello che mi serve.
Ciao dissonance. Volevo rispondere prima ma sono molto occupato in questi giorni. Ho seguito il discorso che avete fatto e vi faccio notare una cosa che mi è saltata subito agli occhi. L'uso delle matrici di permutazione che fate a me sembra sbagliato (dico sembra perchè non vorrei aver capito male io). Nel senso che la moltiplicazione a destra per una matrice di permutazione non scambia le righe ma le colonne.
Nel tuo caso ad ogni passo k moltiplichi la matrice A di quel passo (quella che è stata triangolata fino alla colonna k) a sinistra prima per una Permutazione $Pi$ e poi per una Gauss $M$.
Ora, se consideri:
$A^((k+1))=M^((k))Pi^((k))A^((k))=M^((k))(Pi^((k))A^((k)))$ vuol dire che prendi la matrice del sistema, gli scambi le righe e poi applichi la Gauss.
Se invece (come hai tentato di fare su) consideri:
$A^((k+1))=M^((k))Pi^((k))A^((k))=(M^((k))Pi^((k)))A^((k))$ vuol dire che prendi la matrice del sistema e la moltiplichi direttamente per una matrice che è la Gauss del passo k+1 con le colonne scambiate. E se scambi le colonne di una Gauss (più precisamente se scambi l'indice k con un indice successivo come si fa nella variante del massimo pivot parziale) ottieni ancora una Gauss, in cui semplicemente hai cambiato pivot, però la struttura è la stessa.
Qui non so se ho frainteso io. Però ho l'impressione che tu abbia confuso questo con la variante del massimo pivot totale. In quest'altra infatti si fanno sia scambi di riga che di colonna e allora fare questo giochetto non funziona più per il semplice motivo che anche il vettore delle soluzioni viene permutato.
Un consiglio se, come ho capito, vuoi implementare l'algoritmo (sul pc??). Non è necessario che scrivi nell'algoritmo una procedura che scambi effettivamente le righe della matrice del sistema. Di solito si usano vettori flags, ovvero un vettore di n componenti che inizialmente contiene il numero i alla posizione i. Ad ogni passo invece di scambiare le righe della matrice scambi le componenti del flag e quando esegui materialmente i conti, invece di accedere alla riga i della matrice accedi alla riga flag(i).
Risparmi un bel po' in termini di costo computazionale. Nel caso del pivot totale ti servono 2 vettori flag, uno per le righe e uno per le colonne.
Un altro accorgimento: tutti i calcoli possono essere eseguiti nello stesso spazio di memoria. Man mano che esegui i calcoli modifichi la parte triangolare superiore della matrice iniziale A che alla fine diventerà la U, mentre nella parte inferiore memorizzi man mano che le calcoli le componenti delle matrici di Gauss sovrascrivendo le componenti della matrice iniziale.
Scusa se forse sono stato lapidario, spero nella fretta di non aver scritto scemenze, se non è chiaro qualcosa chiedi pure e appena ho qualche minuto ti rispondo.
Ciao!
Nel tuo caso ad ogni passo k moltiplichi la matrice A di quel passo (quella che è stata triangolata fino alla colonna k) a sinistra prima per una Permutazione $Pi$ e poi per una Gauss $M$.
Ora, se consideri:
$A^((k+1))=M^((k))Pi^((k))A^((k))=M^((k))(Pi^((k))A^((k)))$ vuol dire che prendi la matrice del sistema, gli scambi le righe e poi applichi la Gauss.
Se invece (come hai tentato di fare su) consideri:
$A^((k+1))=M^((k))Pi^((k))A^((k))=(M^((k))Pi^((k)))A^((k))$ vuol dire che prendi la matrice del sistema e la moltiplichi direttamente per una matrice che è la Gauss del passo k+1 con le colonne scambiate. E se scambi le colonne di una Gauss (più precisamente se scambi l'indice k con un indice successivo come si fa nella variante del massimo pivot parziale) ottieni ancora una Gauss, in cui semplicemente hai cambiato pivot, però la struttura è la stessa.
Qui non so se ho frainteso io. Però ho l'impressione che tu abbia confuso questo con la variante del massimo pivot totale. In quest'altra infatti si fanno sia scambi di riga che di colonna e allora fare questo giochetto non funziona più per il semplice motivo che anche il vettore delle soluzioni viene permutato.
Un consiglio se, come ho capito, vuoi implementare l'algoritmo (sul pc??). Non è necessario che scrivi nell'algoritmo una procedura che scambi effettivamente le righe della matrice del sistema. Di solito si usano vettori flags, ovvero un vettore di n componenti che inizialmente contiene il numero i alla posizione i. Ad ogni passo invece di scambiare le righe della matrice scambi le componenti del flag e quando esegui materialmente i conti, invece di accedere alla riga i della matrice accedi alla riga flag(i).
Risparmi un bel po' in termini di costo computazionale. Nel caso del pivot totale ti servono 2 vettori flag, uno per le righe e uno per le colonne.
Un altro accorgimento: tutti i calcoli possono essere eseguiti nello stesso spazio di memoria. Man mano che esegui i calcoli modifichi la parte triangolare superiore della matrice iniziale A che alla fine diventerà la U, mentre nella parte inferiore memorizzi man mano che le calcoli le componenti delle matrici di Gauss sovrascrivendo le componenti della matrice iniziale.
Scusa se forse sono stato lapidario, spero nella fretta di non aver scritto scemenze, se non è chiaro qualcosa chiedi pure e appena ho qualche minuto ti rispondo.
Ciao!
Ps. Non ho sottolineato che da tutto ciò che ho detto segue che effettivamente la fattorizzazione con tecnica di massimo pivot parziale coincide con la fattorizzazione standard della matrice $Pi * A$. Un'osservazione che credo sia bene fare è che in generale la $Pi$ per cui esiste la fattorizzazione di $Pi * A$ non è unica (anzi credo non lo sia mai), quindi quella che emerge dalla tecnica del massimo pivot è solo una possibilità, che numericamente sembra la migliore (totale è meglio di parziale ma costa di più) perchè riesci a contenere gli errori di arrotondamento dovuti a pivot quasi-nulli.
Ciao!!! In effetti questo topic è nato per una cosa e adesso si evolve verso un'altra. Perciò provo a mettere un po' di chiarezza ripartendo da zero.
Quello che mi servirebbe è trovare un algoritmo stabile per la fattorizzazione LU di una matrice A (se applicabile) o di una matrice $Pi$A (A con le righe permutate). Naturalmente nel caso sia necessaria una permutazione delle righe, devo poter essere in grado di conoscerla esattamente al termine dell'algoritmo.
Avevo pensato di usare l'algoritmo di Gauss, che (se applicabile) fattorizza una matrice A in LU. Ora, Gauss "nudo e crudo" non è sempre applicabile e inoltre è poco stabile. Quindi si tratta di applicare una variante del pivot: per semplicità usiamo quella parziale. In sostanza allora, ad ogni passo applichiamo una permutazione delle righe.
Secondo quel brano che ho riportato, alla fine dell'algoritmo avremo calcolato la fattorizzazione LU di una matrice $Pi$A. Domanda: chi è $Pi$? Mi serve saperlo, altrimenti non so cosa ho fatto esattamente. Quale matrice ho fattorizzato?
Se ad esempio dovessi risolvere una famiglia di sistemi, in cui ogni colonna di termini noti dipende dalla soluzione del sistema precedente, conoscendo $Pi$ sarei a cavallo: mi basterebbe risolvere ogni volta $PiA=b$ cioè $LU=Pib$. Ma se non conosco $Pi$, di questa fattorizzazione non me ne faccio nulla.
Quello che mi servirebbe è trovare un algoritmo stabile per la fattorizzazione LU di una matrice A (se applicabile) o di una matrice $Pi$A (A con le righe permutate). Naturalmente nel caso sia necessaria una permutazione delle righe, devo poter essere in grado di conoscerla esattamente al termine dell'algoritmo.
Avevo pensato di usare l'algoritmo di Gauss, che (se applicabile) fattorizza una matrice A in LU. Ora, Gauss "nudo e crudo" non è sempre applicabile e inoltre è poco stabile. Quindi si tratta di applicare una variante del pivot: per semplicità usiamo quella parziale. In sostanza allora, ad ogni passo applichiamo una permutazione delle righe.
Secondo quel brano che ho riportato, alla fine dell'algoritmo avremo calcolato la fattorizzazione LU di una matrice $Pi$A. Domanda: chi è $Pi$? Mi serve saperlo, altrimenti non so cosa ho fatto esattamente. Quale matrice ho fattorizzato?
Se ad esempio dovessi risolvere una famiglia di sistemi, in cui ogni colonna di termini noti dipende dalla soluzione del sistema precedente, conoscendo $Pi$ sarei a cavallo: mi basterebbe risolvere ogni volta $PiA=b$ cioè $LU=Pib$. Ma se non conosco $Pi$, di questa fattorizzazione non me ne faccio nulla.
Sì ok, forse ho infilato troppa roba.
Senti io lo guarderi così: anche se si parla di matrici di permutazione pensale come le permutazioni del gruppo simmetrico in algebra.
Ad passo 1 ad esempio scambi la prima riga con la 4°. Allora hai la permutazione (14). Al secondo passo scambi la seconda riga con la 5°. E allora hai (25). Potrebbero esserci ripetizioni.
Alla fine fai il prodotto di tutte le permutazioni che hai eseguito in ordine inverso. L'idea è che ad ogni passo k non hai una matrice di permutazione qualunque ma una matrice di permutazione che scambia esattamente due righe di cui una è la k-esima e l'altra ha indice maggiore di k.
Il che vuol dire che superato il passo k non toccherai più le righe dalla 1 alla k. Quindi alla fine sai dire esattamente cosa avrai fatto.
Tornando a parlare di matrici succede la stessa cosa, le pui trattare ciascuna come scambi di 2 righe. (se riguardi quello che ti avevo scritto sull'uso di flags era proprio questo quello che intendevo). Alla fine il vettore flag sarà qualcosa del tipo (137264859...) cioè se lo scrivi come prodotto di cicli (come si fa nei corsi di algebra) e lo trasformi in matrice di permutazione ottieni esattamente la $Pi$ che stai cercando. La conclusione è che tratti tutto il procedimento come se fosse un prodotto di cicli (tramite l'isomorfismo che puoi immaginare) e quindi lavori in un ambiente più facile.
Come verificarlo? Basta pensare a quello che succede al vettore dei termini noti, in cui le cose sono più facili visto che lo scambio è tra componenti è ogni moltiplicazione per una Gauss risulta in (n-k) somme tra le componenti del vettore.
Non so se sono stato abbastanza chiaro.
Credo fermamente che tentare di estrarre la matrice $Pi$ dal prodotto $M^nPi^n....M^1Pi^1$ sia (nel caso generale) un'impresa impossibile o comunque un po' troppo lunga.
Senti io lo guarderi così: anche se si parla di matrici di permutazione pensale come le permutazioni del gruppo simmetrico in algebra.
Ad passo 1 ad esempio scambi la prima riga con la 4°. Allora hai la permutazione (14). Al secondo passo scambi la seconda riga con la 5°. E allora hai (25). Potrebbero esserci ripetizioni.
Alla fine fai il prodotto di tutte le permutazioni che hai eseguito in ordine inverso. L'idea è che ad ogni passo k non hai una matrice di permutazione qualunque ma una matrice di permutazione che scambia esattamente due righe di cui una è la k-esima e l'altra ha indice maggiore di k.
Il che vuol dire che superato il passo k non toccherai più le righe dalla 1 alla k. Quindi alla fine sai dire esattamente cosa avrai fatto.
Tornando a parlare di matrici succede la stessa cosa, le pui trattare ciascuna come scambi di 2 righe. (se riguardi quello che ti avevo scritto sull'uso di flags era proprio questo quello che intendevo). Alla fine il vettore flag sarà qualcosa del tipo (137264859...) cioè se lo scrivi come prodotto di cicli (come si fa nei corsi di algebra) e lo trasformi in matrice di permutazione ottieni esattamente la $Pi$ che stai cercando. La conclusione è che tratti tutto il procedimento come se fosse un prodotto di cicli (tramite l'isomorfismo che puoi immaginare) e quindi lavori in un ambiente più facile.
Come verificarlo? Basta pensare a quello che succede al vettore dei termini noti, in cui le cose sono più facili visto che lo scambio è tra componenti è ogni moltiplicazione per una Gauss risulta in (n-k) somme tra le componenti del vettore.
Non so se sono stato abbastanza chiaro.
Credo fermamente che tentare di estrarre la matrice $Pi$ dal prodotto $M^nPi^n....M^1Pi^1$ sia (nel caso generale) un'impresa impossibile o comunque un po' troppo lunga.
Si, il fatto che ci sia isomorfismo tra matrici di permutazione e permutazioni lo sapevo.Vediamo se ho capito quello che dici tu, a parole terra-terra: "che te ne frega di ricavare esattamente la permutazione dal prodotto $M_nPi_n...$, ti memorizzi gli scambi che fai a ogni passo e sei in grado di ricostruire tutto quello che hai fatto". Giusto? Inoltre dici che secondo te usare Gauss per la fattorizzazione LU a meno di matrici di permutazione non va bene, sempre perché:
Ho capito bene? Se sì, su questo devo ragionarci un po' (stasera purtroppo sono troppo cotto però).
Nel frattempo pongo una domanda a cui penso sarà più facile rispondere: che algoritmo numericamente stabile posso usare per la fattorizzazione LU?
Credo fermamente che tentare di estrarre la matrice Π dal prodotto MnΠn....M1Π1 sia (nel caso generale) un'impresa impossibile o comunque un po' troppo lunga.
Ho capito bene? Se sì, su questo devo ragionarci un po' (stasera purtroppo sono troppo cotto però).
Nel frattempo pongo una domanda a cui penso sarà più facile rispondere: che algoritmo numericamente stabile posso usare per la fattorizzazione LU?
Nono, forse mi sono espresso male. Il metodo di Gauss, che io sappia, è il migliore per la fatt. LU nel caso di matrici generiche. Devi necessariamente prevedere un sistema di pivot perlomeno parziali. Esistono le cosidette tecniche compatte di cui so poco, ma non penso che ci ricavi più di tanto. In ogni caso l'algoritmo di risoluzione di un sistema lineare con Gauss in generale non è una buona scelta numerica, si studia più che altro per motivi storici. Ci sono particolari matrici per cui la fattorizzazione di Gauss conviene, ma in generale si usano altre tecniche. Questo stando alle informazioni che ho io.
Volendo implementare l'algoritmo di Gauss c'è da scegliere tra pivoting parziale e totale. Il secondo assicura una maggiorazione dell'errore molto inferiore al primo. Tieni presente che in entrambi i casi la funzione che maggiora ha andamento esponenziale. C'è una formula che ora non ricordo, con una radice e una produttoria. Credo la trovi in ogni manuale di analisi numerica. Per curiosità potresti anche dare un'occhiata alla congettura di Wilkinson.
Come scegliere tra parziale e totale? Prova a pensare ad un modo per "stimare" la necessità e il costo dell'uno o dell'altro leggendo i valori delle entrate della matrice. Deve esserti ben chiaro però che non esiste il metodo perfetto, ciò che guadagni in sbabilità lo perdi in costo computazionale e/o di memorizzazione. Inoltre c'è il malcondizionamento sempre pronto dietro l'angolo e stimare quello non è semplice dal punto di vista computazionale nel caso di matrici generiche.
Volendo implementare l'algoritmo di Gauss c'è da scegliere tra pivoting parziale e totale. Il secondo assicura una maggiorazione dell'errore molto inferiore al primo. Tieni presente che in entrambi i casi la funzione che maggiora ha andamento esponenziale. C'è una formula che ora non ricordo, con una radice e una produttoria. Credo la trovi in ogni manuale di analisi numerica. Per curiosità potresti anche dare un'occhiata alla congettura di Wilkinson.
Come scegliere tra parziale e totale? Prova a pensare ad un modo per "stimare" la necessità e il costo dell'uno o dell'altro leggendo i valori delle entrate della matrice. Deve esserti ben chiaro però che non esiste il metodo perfetto, ciò che guadagni in sbabilità lo perdi in costo computazionale e/o di memorizzazione. Inoltre c'è il malcondizionamento sempre pronto dietro l'angolo e stimare quello non è semplice dal punto di vista computazionale nel caso di matrici generiche.
No...In realtà chi si è espresso male sono io. Ho fatto trecento domande a raffica, una sovrapposta all'altra e non ho fatto capire niente! In questo momento, non ho il problema di scegliere tra pivoting parziale e totale, e nemmeno tra algoritmo di Gauss o altro.
Il mio problema è che non capisco come l'algoritmo di Gauss con uno dei due pivoting possa restituirci la fattorizzazione LU di una matrice. Infatti l'algoritmo senza pivoting, se non si inceppa(quindi se i minori principali di testa sono non nulli), ci fa uno scherzo del genere: al termine dell'algoritmo abbiamo $E_(n-1)...E_(1)A=U$. Invertendo la trafila di matrici elementari otteniamo una matrice $L$, da cui $A=LU$.
Questo era senza nessun pivoting.
Se invece introduciamo una strategia di pivoting, inframmezzate tra le matrici elementari $E_(n-1)...E_(1)$ avremo delle matrici di permutazione (che poi saranno tutte scambi semplici, pacifico). Ma allora da dove la ricaviamo la matrice $L$?
Il mio problema è che non capisco come l'algoritmo di Gauss con uno dei due pivoting possa restituirci la fattorizzazione LU di una matrice. Infatti l'algoritmo senza pivoting, se non si inceppa(quindi se i minori principali di testa sono non nulli), ci fa uno scherzo del genere: al termine dell'algoritmo abbiamo $E_(n-1)...E_(1)A=U$. Invertendo la trafila di matrici elementari otteniamo una matrice $L$, da cui $A=LU$.
Questo era senza nessun pivoting.
Se invece introduciamo una strategia di pivoting, inframmezzate tra le matrici elementari $E_(n-1)...E_(1)$ avremo delle matrici di permutazione (che poi saranno tutte scambi semplici, pacifico). Ma allora da dove la ricaviamo la matrice $L$?
Accedi a tutti gli appunti
Tutor AI: studia meglio e in meno tempo