[Ottimizzazione] Risoluzione sistema lineare

[HeroGian]

Salve, devo implementare per esercizio un programmino matlab che risolva un sistema lineare tramite minimizzazione di una funzione quadratica con un metodo del gradiente con:

-direzione calcolata tramite Steepest descent
-steplength calcolato dalla regola di armijo

function y = f(A, x, b)
	y = (1/2)*x'*A*x -b'*x;
end

function alfa = armijo(A, x, b, g, s, beta, sigma)
	alfa = s;

	while f(A, x -alfa*g, b) > f(A, x, b)+sigma*alfa*g'*g
		alfa = beta*alfa;
	end
end

function x = steepest(A, b, max_it, toll)
	n = length(A);
	x_prec = rand(n, 1);
	g = A*x_prec-b;

	for k=1 : max_it
		if norm(g) < toll*norm(b)
		    disp('condizione sulla norma raggiunta');
		    break;
		end
		alfa = armijo(A, x_prec, b, g, 1, 0.5, 1e-4);
		x = x_prec -alfa*g;
		g = A*x-b;
		x_prec = x;
	end
end

mi pare uguale a com'è dal punto di vista teorico: tramite armijo mi calcolo lo steplength e costruisco la successione di vettori in direzione opposta al gradiente.. solo che dopo un po di iterazioni sembra che le componenti di x vadano in overflow..

[HeroGian]

Problema risolto, in realtà l'algoritmo funziona.. ma non avevo notato che le matrici con cui stavo testando non erano simmetriche definite positive