[Algoritmo MATLAB] Metodo del Gradiente (rapida scesa)
Salve a tutti. Un docente mi ha affidato questo esercizio:
"Assegnato il sistema lineare $Ax = b$ , $ A \in \mathbb R^{n \times n} $ simmetrica e definita positiva,
scrivere un sottoprogramma Matlab per il calcolo della sua soluzione usando il
metodo del gradiente.
Per gli associati problemi di minimizzazione unidimensionale utilizzare una
ricerca su un intervallo discretizzato.
Verificare la correttezza del sottoprogramma mediante opportuni problemi test"
Io ho fatto il seguente algoritmo:
Ma i risultati sballano, ci deve essere qualche errore! Mi aiutate???
"Assegnato il sistema lineare $Ax = b$ , $ A \in \mathbb R^{n \times n} $ simmetrica e definita positiva,
scrivere un sottoprogramma Matlab per il calcolo della sua soluzione usando il
metodo del gradiente.
Per gli associati problemi di minimizzazione unidimensionale utilizzare una
ricerca su un intervallo discretizzato.
Verificare la correttezza del sottoprogramma mediante opportuni problemi test"
Io ho fatto il seguente algoritmo:
function[x,iter,err]=gradedue(A,b,x0,h,N,nmax,toll)
r=b-A*x0;
s=linspace(-h,h,N);
f=zeros(size(s));
iter=0;
r0=norm(r);
err=norm(r);
x=x0;
while err>toll&&iter<nmax
iter=iter+1;
for i=1:N
f(i)=(0.5*((x+s(i)*r)'*A*(x+s(i)*r)))-((x+s(i)*r)'*b);
end
k=min(f);
x=x+k*r;
r=b-A*x;
err=norm(r)/r0;
end
return
Ma i risultati sballano, ci deve essere qualche errore! Mi aiutate???
Risposte
Per prima cosa, utilizza il tag code per rendere il codice un po' più leggibile.
Veniamo quindi al tuo problema. Devo studiarmi il metodo un po' più nel dettaglio per dirti dove sbagli, ma una cosa molto utile quando si cerca di verificare la correttezza del metodo del gradiente è quella di fare un grafico dell'andamento dell'errore in funzione delle iterazioni. Com'è questo grafico? L'errore decresce? Con quale "velocità"? E' infatti anche possibile che l'errore sia nella scelta di h, toll o nmax...
L'unica cosa che per ora osservo è che scegli il parametro del problema di minimizzazione unidimensionale in \( [-h, h] \). A prima vista non mi sembra però sensato scegliere un parametro negativo.. Sarebbe come decidere che la scelta migliore è quella di allontanarti dal risultato corretto. Sbaglio?
Veniamo quindi al tuo problema. Devo studiarmi il metodo un po' più nel dettaglio per dirti dove sbagli, ma una cosa molto utile quando si cerca di verificare la correttezza del metodo del gradiente è quella di fare un grafico dell'andamento dell'errore in funzione delle iterazioni. Com'è questo grafico? L'errore decresce? Con quale "velocità"? E' infatti anche possibile che l'errore sia nella scelta di h, toll o nmax...
L'unica cosa che per ora osservo è che scegli il parametro del problema di minimizzazione unidimensionale in \( [-h, h] \). A prima vista non mi sembra però sensato scegliere un parametro negativo.. Sarebbe come decidere che la scelta migliore è quella di allontanarti dal risultato corretto. Sbaglio?
Scusami per non aver usato il tag quote, dopo modifico il messaggio.
Per il -h...alla fine io prendo un vettore e metto le sue componenti nel formulone che mi da f, dalla cui componente minima (k) mi ricavo la "posizione" x...che dovrebbe essere la più vicina alla soluzione. Perchè non dovrebbero esserci componenti negative nel vettore? Il minimo assoluto (per meglio dire l'approssimazione del minimo assoluto) che cerco, non potrebbe essere anche negativa?
Per il -h...alla fine io prendo un vettore e metto le sue componenti nel formulone che mi da f, dalla cui componente minima (k) mi ricavo la "posizione" x...che dovrebbe essere la più vicina alla soluzione. Perchè non dovrebbero esserci componenti negative nel vettore? Il minimo assoluto (per meglio dire l'approssimazione del minimo assoluto) che cerco, non potrebbe essere anche negativa?
Non è il tag quote da usare, ma il tag code. Ho comunque modificato il tuo post in modo che usasse il tag corretto.
Io ho sempre visto considerare solo parametri positivi. Dopotutto stai decidendo di quanto muoverti lungo una particolare direzione e a parametri negativi corrispondono movimenti in direzione opposta. Hai provato a vedere se l'errore decresce abbastanza velocemente? Hai verificato la formula da minimizzare? Puoi postare gli input che hai dato alla funzione e quelli che ritieni debbano essere gli output?
Io ho sempre visto considerare solo parametri positivi. Dopotutto stai decidendo di quanto muoverti lungo una particolare direzione e a parametri negativi corrispondono movimenti in direzione opposta. Hai provato a vedere se l'errore decresce abbastanza velocemente? Hai verificato la formula da minimizzare? Puoi postare gli input che hai dato alla funzione e quelli che ritieni debbano essere gli output?
Una prova che ho fatto è con il seguente esempio:
Mentre utilizzando
Oppure:
e ancora:
Ho anche provato a cambiare tipo di algoritmo, cioè creare un sottoprogramma che mi calcola min(f) per poi richiamarlo in un algoritmo generale, ma il risultato non cambia.
Come posso vedere se l'errore decresce abbastanza velocemente? La formula da minimizzare mi sembra corretta. Magari per il vettore ad estremi positivi hai ragione tu.
A=[2,-1,0,0;-1,2,-1,0;0,-1,2,-1;0,0,-1,-2]
x0=[1;1;1;1]
b=[2;3;4;4;]
[x,iter,err]=gradedue(A,b,x0,100,11,3,0.003)
Risultato:
x =
1.0e+080 *
0.6325
-1.0541
-2.7407
-8.0112
iter =
3
err =
2.2219e+080Mentre utilizzando
[x,iter,err]=gradedue(A,b,x0,100,11,100,0.00003)
x =
NaN
NaN
NaN
NaN
iter =
5
err =
NaNOppure:
[x,iter,err]=gradedue(A,b,x0,100,11,12,0.003)
x =
0
0
0
iter =
12
err =
1
e ancora:
[x,iter,err]=gradedue(A,b,x0,100,11,1200,0.003)
x =
0
0
0
iter =
1200
err =
1
Ho anche provato a cambiare tipo di algoritmo, cioè creare un sottoprogramma che mi calcola min(f) per poi richiamarlo in un algoritmo generale, ma il risultato non cambia.
Come posso vedere se l'errore decresce abbastanza velocemente? La formula da minimizzare mi sembra corretta. Magari per il vettore ad estremi positivi hai ragione tu.
Farei un piccolo ripasso della teoria. Per prima cosa l'obiettivo della tua funzione è quello di trovare la soluzione al sistema lineare \( A\,x = b \) minimizzando la funzione \( f(x) = || b - A\,x ||^2. \) Per trovare questo minimo parti da una approssimazione e ti muovi nella direzione di massima discesa nell'intorno di tale punto, cioè nella direzione opposta al gradiente. Se lo spostamento è abbastanza piccolo, la funzione dovrebbe decrescere. Aggiorni quindi la tua approssimazione \( x_n \) con la formula \( x_{n+1} = x_n - \lambda\,\nabla\,f\,(x_n) \) per un qualche \( \lambda > 0. \) E' corretto fino a qui? La tua teoria è diversa?
A questo punto dobbiamo trovare il gradiente di
\[ || b - A\,x ||^2 = \langle b - A\,x, b - A\,x \rangle = \langle b, b \rangle - 2\,\langle A\,x, b \rangle + \langle A\,x, A\,x \rangle. \]
Se non ho sbagliato i calcoli dovrebbe essere
\[ \nabla\,f\,(x) = - 2\,A^T\,b + 2\,(A^T\,A)\,x = 2\,A^T\,( A\,x - b ). \]
Questa formula non mi sembra uguale a quella che hai utilizzato nel tuo codice. Stai quindi facendo qualcosa di diverso?
Non capisco poi il seguente codice:
Perché consideri l'errore relativo a r0 che è costante dalla prima iterazione?
Perché passi gli argomenti \(h = 100\) e \( N = 11 \) alla funzione? In generale lo spostamento lungo il gradiente negativo deve essere piccolo. Nel tuo codice consideri valori per lo spostamento molto alti.. Perché hai scelto spesso un numero di iterazioni così basso?
A questo punto dobbiamo trovare il gradiente di
\[ || b - A\,x ||^2 = \langle b - A\,x, b - A\,x \rangle = \langle b, b \rangle - 2\,\langle A\,x, b \rangle + \langle A\,x, A\,x \rangle. \]
Se non ho sbagliato i calcoli dovrebbe essere
\[ \nabla\,f\,(x) = - 2\,A^T\,b + 2\,(A^T\,A)\,x = 2\,A^T\,( A\,x - b ). \]
Questa formula non mi sembra uguale a quella che hai utilizzato nel tuo codice. Stai quindi facendo qualcosa di diverso?
Non capisco poi il seguente codice:
err=norm(r)/r0;
Perché consideri l'errore relativo a r0 che è costante dalla prima iterazione?
Perché passi gli argomenti \(h = 100\) e \( N = 11 \) alla funzione? In generale lo spostamento lungo il gradiente negativo deve essere piccolo. Nel tuo codice consideri valori per lo spostamento molto alti.. Perché hai scelto spesso un numero di iterazioni così basso?
Dunque!
Ti spiego cosa dice la "mia teoria":
come dici tu, questo è un problema di minimizzazione, cerco le soluzione del sistema lineare Ax=b con A simmetrica e definita positiva come i punti di minimo della funzione:
$ f(x)=1/2((x)^(T)Ax)-(x^T)b $
La derivata sarà Ax-b.
Considero $r^(k)=b-Ax^(k)$ dove k è il k-appesimo passaggio. Se $r^k$ è diverso da 0 allora esiste $v^k$ soluzione del problema di minimizzazione (nel mio caso quello del vettore da cui prendo le componenti) tale che:$f(x^(k)+v^(k)r^k)<=f(x^(k)+vr^k)$ per ogni v.
Il metodo del gradiente mi dice che $x^(k+1)=x^(k)+v^(k)r^(k)$
Questo per la teoria.
Per quanto riguarda quel codice mi è stato consigliato di utilizzarlo così. Ho provato con altri programmi e non mi da errori.
Per quanto riguarda le domande sugli esempi ho preso a caso quei dati per testare il metodo. Se funzionasse bene dovrebbe andare anche con quei dati no?
Ti spiego cosa dice la "mia teoria":
come dici tu, questo è un problema di minimizzazione, cerco le soluzione del sistema lineare Ax=b con A simmetrica e definita positiva come i punti di minimo della funzione:
$ f(x)=1/2((x)^(T)Ax)-(x^T)b $
La derivata sarà Ax-b.
Considero $r^(k)=b-Ax^(k)$ dove k è il k-appesimo passaggio. Se $r^k$ è diverso da 0 allora esiste $v^k$ soluzione del problema di minimizzazione (nel mio caso quello del vettore da cui prendo le componenti) tale che:$f(x^(k)+v^(k)r^k)<=f(x^(k)+vr^k)$ per ogni v.
Il metodo del gradiente mi dice che $x^(k+1)=x^(k)+v^(k)r^(k)$
Questo per la teoria.
Per quanto riguarda quel codice mi è stato consigliato di utilizzarlo così. Ho provato con altri programmi e non mi da errori.
Per quanto riguarda le domande sugli esempi ho preso a caso quei dati per testare il metodo. Se funzionasse bene dovrebbe andare anche con quei dati no?
A quanto pare è stato utile ridefinire la teoria in modo da capire che stavamo pensando a cose diverse. Il metodo del gradiente ha in effetti un'utilità molto più generale di questi due casi particolari.
I casi particolare da testare non vanno scelti proprio così a caso, e questo vale ancora meno per i parametri come h, N, mmax e toll. Per prima cosa il problema verificato deve rispettare le condizioni del tuo problema originario. \(A\) deve quindi prima di tutto essere simmetrica e definita positiva. La matrice \(A\) del tuo esempio non è definita positiva. Ha infatti come autovalori (calcolati con Octave): \(-2.24995\), \(0.66563\), \(2.12046\) e \(3.46386\). Probabilmente hai sbagliato a scrivere il segno del \(-2\) nell'ultima riga. Se ci fosse infatti stato un \(2\) al posto del \(-2\) avresti avuto una matrice definita positiva.
Dopodiché è utile, MOLTO utile, conoscere il risultato del problema che stai verificando. Ti permette infatti di verificare che passando il risultato corretto come valore iniziale, la funzione restituisce questo risultato corretto (o almeno una sua approssimazione molto vicina). A questo proposito potrebbe convenire partire da problemi molto semplici come \(A = I\) in modo da verificare che il problema funzioni almeno in casi semplici come questo.
Il tuo programma è poi diviso in due parti: il metodo del gradiente e la ricerca lineare del parametro migliore per avanzare con tale metodo. Conviene a mio parere testare le due parti separatamente in modo da essere certo che il programma stia funzionando correttamente in ogni sua parte. Puoi per prima cosa definire il metodo del gradiente in modo che riceva un coefficiente di avanzamento fisso (che non faccia cioè la ricerca lineare per ottimizzare questo coefficiente). Conviene sceglierlo abbastanza basso e verificare (plottando l'evoluzione dell'errore) che l'errore decresca. Alzando questo coefficiente l'errore dovrebbe scendere più velocemente ma c'è il rischio che l'errore possa divergere o che oscilli. Quando questo valore è molto basso la convergenza c'è, ma è potenzialmente molto lenta. In molte situazione pratiche la ricerca binaria non si fa. Al suo posto si sceglie un valore iniziale per questo coefficiente e si esegue il metodo. Se questo non converge, allora il coefficiente si diminuisce e la ricerca del risultato riparte. Ti consiglio di provare diversi parametri per vedere l'effetto che questa scelta ha sul comportamento del metodo.
Una volta che hai verificato che il tuo metodo funziona senza l'ottimizzazione del parametro, prova ad inserire tale ricerca lineare. Cerca di usare quello che hai imparato nella fase precedente per scegliere su quali parametri effettuare la ricerca.
Riguardo alla tolleranza e al numero massimo di iterazioni invece, usa numeri abbastanza alti per iniziare. Questo metodo converge con grande lentezza..
I casi particolare da testare non vanno scelti proprio così a caso, e questo vale ancora meno per i parametri come h, N, mmax e toll. Per prima cosa il problema verificato deve rispettare le condizioni del tuo problema originario. \(A\) deve quindi prima di tutto essere simmetrica e definita positiva. La matrice \(A\) del tuo esempio non è definita positiva. Ha infatti come autovalori (calcolati con Octave): \(-2.24995\), \(0.66563\), \(2.12046\) e \(3.46386\). Probabilmente hai sbagliato a scrivere il segno del \(-2\) nell'ultima riga. Se ci fosse infatti stato un \(2\) al posto del \(-2\) avresti avuto una matrice definita positiva.
Dopodiché è utile, MOLTO utile, conoscere il risultato del problema che stai verificando. Ti permette infatti di verificare che passando il risultato corretto come valore iniziale, la funzione restituisce questo risultato corretto (o almeno una sua approssimazione molto vicina). A questo proposito potrebbe convenire partire da problemi molto semplici come \(A = I\) in modo da verificare che il problema funzioni almeno in casi semplici come questo.
Il tuo programma è poi diviso in due parti: il metodo del gradiente e la ricerca lineare del parametro migliore per avanzare con tale metodo. Conviene a mio parere testare le due parti separatamente in modo da essere certo che il programma stia funzionando correttamente in ogni sua parte. Puoi per prima cosa definire il metodo del gradiente in modo che riceva un coefficiente di avanzamento fisso (che non faccia cioè la ricerca lineare per ottimizzare questo coefficiente). Conviene sceglierlo abbastanza basso e verificare (plottando l'evoluzione dell'errore) che l'errore decresca. Alzando questo coefficiente l'errore dovrebbe scendere più velocemente ma c'è il rischio che l'errore possa divergere o che oscilli. Quando questo valore è molto basso la convergenza c'è, ma è potenzialmente molto lenta. In molte situazione pratiche la ricerca binaria non si fa. Al suo posto si sceglie un valore iniziale per questo coefficiente e si esegue il metodo. Se questo non converge, allora il coefficiente si diminuisce e la ricerca del risultato riparte. Ti consiglio di provare diversi parametri per vedere l'effetto che questa scelta ha sul comportamento del metodo.
Una volta che hai verificato che il tuo metodo funziona senza l'ottimizzazione del parametro, prova ad inserire tale ricerca lineare. Cerca di usare quello che hai imparato nella fase precedente per scegliere su quali parametri effettuare la ricerca.
Riguardo alla tolleranza e al numero massimo di iterazioni invece, usa numeri abbastanza alti per iniziare. Questo metodo converge con grande lentezza..
Allora!
Evidentemente ho sbagliato a scrivere si! Per i 2 programmi in cui è diviso il programma li ho già provati separatamente e funzionano. Perchè in origine avevo fatto il metodo del gradiente diretto e andava. Poi ho aggiunto l'altro, provandolo prima da solo e andava.
Ho chiesto lumi riguardo alla perplessità che sollevavi tu, sulla scelta del -h. Mi è stato detto che è vero e poi mi è stato detto (cosa che non ho capito però) che per il calcolo dell'iterata x^{k+1} conviene comunque prendere dei valori negativi [-k,h] con k
Questi esempi possono andare?
Le soluzioni sono in entrambi i casi un vettore unitario.
Ah e comunque a questo punto come modifico la parte s=linspace(-h,h,N) ?
Poi senti, io questo algoritmo lo dovrei finire entro il 14. Avresti tempo per fare una sessioncina skype-msn?
Evidentemente ho sbagliato a scrivere si! Per i 2 programmi in cui è diviso il programma li ho già provati separatamente e funzionano. Perchè in origine avevo fatto il metodo del gradiente diretto e andava. Poi ho aggiunto l'altro, provandolo prima da solo e andava.
Ho chiesto lumi riguardo alla perplessità che sollevavi tu, sulla scelta del -h. Mi è stato detto che è vero e poi mi è stato detto (cosa che non ho capito però) che per il calcolo dell'iterata x^{k+1} conviene comunque prendere dei valori negativi [-k,h] con k
Questi esempi possono andare?
A=[1,2,0;2,1,2;0,2,1]
A =
1 2 0
2 1 2
0 2 1
b=[3;5;3]
b =
3
5
3
x0=[0;0;0]
x0 =
0 0 0A = [ 4 -1 0 0 -1 0 0 0 0 0 0 0 ;
-1 4 -1 0 0 -1 0 0 0 0 0 0 ;
0 -1 4 -1 0 0 -1 0 0 0 0 0 ;
0 0 -1 4 0 0 0 -1 0 0 0 0 ;
-1 0 0 0 4 -1 0 0 -1 0 0 0 ;
0 -1 0 0 -1 4 -1 0 0 -1 0 0 ;
0 0 -1 0 0 -1 4 -1 0 0 -1 0 ;
0 0 0 -1 0 0 -1 4 -1 0 0 -1 ;
0 0 0 0 -1 0 0 0 4 -1 0 0 ;
0 0 0 0 0 -1 0 0 -1 4 -1 0 ;
0 0 0 0 0 0 -1 0 0 -1 4 -1 ;
0 0 0 0 0 0 0 -1 0 0 -1 4 ]
b = [ 2 1 1 2 1 0 0 1 2 1 1 2 ]
x0=zeros(12,1)Le soluzioni sono in entrambi i casi un vettore unitario.
Ah e comunque a questo punto come modifico la parte s=linspace(-h,h,N) ?
Poi senti, io questo algoritmo lo dovrei finire entro il 14. Avresti tempo per fare una sessioncina skype-msn?
Secondo me ci sei quasi, e riusciamo a metterlo a posto prima del 14. In questi giorni sono molto occupato ma se non abbiamo ancora trovato il problema entro la fine della settimana cerco di trovare il momento per sentirci su skype o msn.
Mi sembra che i due casi possano andare bene. La ricerca la puoi modificare in molti modi. Se te ne freghi della questione degli errori di arrotondamento puoi effettuare la ricerca in [0, h]. Altrimenti puoi scegliere un valore abbastanza basso (la tolleranza ad esempio) e usarlo per aggiustare il risultato in negativo cercando ad esempio in [-toll, h].
Mi sembra che i due casi possano andare bene. La ricerca la puoi modificare in molti modi. Se te ne freghi della questione degli errori di arrotondamento puoi effettuare la ricerca in [0, h]. Altrimenti puoi scegliere un valore abbastanza basso (la tolleranza ad esempio) e usarlo per aggiustare il risultato in negativo cercando ad esempio in [-toll, h].
Gentilissimo, speriamo che tu abbia ragione!
Sto facendo test con il primo esempio del mio ultimo post. Ho modificato l'algoritmo mettendo linspace[1,h,N]...
E il risultato non cambia se metto -toll al posto di 1 nell' linspace...
Sto facendo test con il primo esempio del mio ultimo post. Ho modificato l'algoritmo mettendo linspace[1,h,N]...
A=[1,2,0;2,1,2;0,2,1]
A =
1 2 0
2 1 2
0 2 1
>> b=[3;5;3]
b =
3
5
3
>> x0=[0;0;0]
x0 =
0
0
0
>> gradedueprov(A,b,x0,50,5,50,0.005)
ans =
-Inf
-Inf
-InfE il risultato non cambia se metto -toll al posto di 1 nell' linspace...
Metti valori troppo alti di \(h\) e troppo bassi di \(N\)... In generale lo spostamento lungo il gradiente negativo deve essere piccolo per non far divergere il metodo. Sarebbe molto meglio avessi scelto \( h = 5 \) e \( N = 50 \)...
Scusatemi il doppio post, ma porto buone nuove.
Allora, mi sono accorto di un errore nell'algoritmo! In quello che c'è sul qui sul forum, io metto k=min(f) che non va bene!
Ecco come ho modificato:
Adesso provo e cedo i risultati.
Ecco! Direi che va meglio...anche se c'è sempre qualche cosa da affinare...
Come posso avvicinarmi al risultato esatto?
Allora, mi sono accorto di un errore nell'algoritmo! In quello che c'è sul qui sul forum, io metto k=min(f) che non va bene!
Ecco come ho modificato:
function[x,iter,err]=gradedueprov(A,b,x0,h,N,nmax,toll)
r=b-A*x0;
s=linspace(-toll,h,N);
f=zeros(size(s));
iter=0;
r0=norm(r);
err=norm(r);
x=x0;
while err>toll&&iter<nmax
iter=iter+1;
for i=1:N
f(i)=(0.5*((x+s(i)*r)'*A*(x+s(i)*r)))-((x+s(i)*r)'*b);
end
[Y,I]=min(f);
x=x+s(I)*r;
r=b-A*x;
err=norm(r)/r0;
end
returnAdesso provo e cedo i risultati.
Ecco! Direi che va meglio...anche se c'è sempre qualche cosa da affinare...
x,iter,err]=gradedueprov(A,b,x0,5,100,100,0.0001)
x =
1.0e+046 *
-0.8732
3.7381
-0.8732
iter =
100
err =
1.4246e+046Come posso avvicinarmi al risultato esatto?
Accedi a tutti gli appunti
Tutor AI: studia meglio e in meno tempo