solve_nonl_CG_num.m
33bba04159d6master

solve_nonl_CG_num.m
View Options

			% x is a vector
			% f = f(x)
			% grad = f(i, x) = df / dx_i
			% hess = f(i, n, x) = d^2 f/(dx_i * dx_k)
			% N = size of X
			function X = solve_nonl_CG_num(f, N, hint)

			if length(hint) == 0
			X = rand(N,1);
			X = X / norm(X);
			else
			X = hint;
			end

			r = -grad(f, X);
			d = r;
			limit = 100;
			alimit = 100;
			k = 0;
			while norm(r) > eps && k < limit
			alpha = 2 * eps / norm(d); % sure \|alpha * d\| > eps
			ak = 0;
			while norm(alpha * d) > eps && ak < alimit
			alpha = - (grad(f, X)' * d) / (d' * hess(f, X) * d);
			X = X + alpha * d;
			ak = ak + 1;
			%fprintf("\|alpha * d\|: %g \n", norm(alpha * d))
			end
			r_old = r(:); % copy r
			r = -grad(f, X);
			beta = (r' * r) / (r_old' * r_old);
			d = r + beta * d;

			k = k + 1;
			fprintf("Norm r: %g\n", norm(r))
			end

			if k >= limit
			warning(sprintf("Solve Non-linear limit reached: \|r\| = %g", norm(r)))
			end
			end

			function g = grad(f, X)
			dx = 0.000001;
			g = zeros(length(X),1);
			for i = 1:length(X)
			Xp = X(:);
			Xp(i) = Xp(i) + dx / 2;
			Xm = X(:);
			Xm(i) = Xm(i) - dx / 2;
			g(i) = (f(Xp) - f(Xm)) / dx;
			end
			end

			function h = hess(f, X)
			dx = 0.000001;
			h = zeros(length(X));
			for i = 1:length(X)
			h(i,:) = grad(f,X)';
			end
			end