Algoritmo di Gauss e fattorizzazione LU

[Megan00b]

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.

[dissonance]

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.

[Megan00b]

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.

[dissonance]

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é:

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?

[Megan00b]

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.

[dissonance]

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$?