steepest descent


architect
04-16-2006, 08:52 AM
Hello All,

My name is Kostas and I am new in the forum. I am new in Matlab as well. I am trying to get some advice in executing a code that I found on a book "Numerical Methods using Matlab". The code is as follows:

function [xo,fo] = opt_steep(f,x0,TolX,TolFun,alpha0,MaxIter)
% minimize the function f by the steepest descent method.
%input: f = ftn to be given as a string ’f’
% x0 = the initial guess of the solution
%output: x0 = the minimum point reached
% f0 = f(x(0))
if nargin < 6, MaxIter = 100; end %maximum # of iteration
if nargin < 5, alpha0 = 10; end %initial step size
if nargin < 4, TolFun = 1e-8; end %|f(x)| < TolFun wanted
if nargin < 3, TolX = 1e-6; end %|x(k)- x(k - 1)|<TolX wanted
x = x0; fx0 = feval(f,x0); fx = fx0;
alpha = alpha0; kmax1 = 25;
warning = 0; %the # of vain wanderings to find the optimum step size
for k = 1: MaxIter
g = grad(f,x); g = g/norm(g); %gradient as a row vector
alpha = alpha*2; %for trial move in negative gradient direction
fx1 = feval(f,x - alpha*2*g);
for k1 = 1:kmax1 %find the optimum step size(alpha) by line search
fx2 = fx1; fx1 = feval(f,x-alpha*g);
if fx0 > fx1+TolFun & fx1 < fx2 - TolFun %fx0 > fx1 < fx2
den = 4*fx1 - 2*fx0 - 2*fx2; num = den - fx0 + fx2; %Eq.(7.1.5)
alpha = alpha*num/den;
x = x - alpha*g; fx = feval(f,x); %Eq.(7.1.9)
break;
else alpha = alpha/2;
end
end
if k1 >= kmax1, warning = warning + 1; %failed to find optimum step size
else warning = 0;
end
if warning >= 2|(norm(x - x0) < TolX&abs(fx - fx0) < TolFun), break; end
x0 = x; fx0 = fx;
end
xo = x; fo = fx;
if k == MaxIter, fprintf(’Just best in %d iterations’,MaxIter), end

%nm714
f713 = inline('x(1)*(x(1) - 4 - x(2)) + x(2)*(x(2)- 1)','x');
x0 = [0 0], TolX = 1e-4; TolFun = 1e-9; alpha0 = 1; MaxIter = 100;
[xo,fo] = opt_steep(f713,x0,TolX,TolFun,alpha0,MaxIter)

The problem occurs on line 15 where the grad is not defined as the error message says in matlab. I would appreciate your help, if somebody more familiar can spot out where the problem is.

Many Thanks

BRs

Kostas

VadimW
12-28-2007, 02:56 PM
%calculate gradient of the function
function gradf=grad(f,x)
lengthx=length(x); % size of input vector
gradf=zeros(size(x)); % intialize the result
for k=1:lengthx % run all over the parameters and calculate the derivative
if x(k)==0 dx=0.001;
else dx=abs(x(k)/100); % take the step size as 1% of the value
end;
diffvect=zeros(size(x)); %zero all the variables in the addition vector except at k
diffvect(k)=dx;
% calculate derivative at k direction
gradf(k)=(feval(f,x+diffvect)-feval(f,x))/dx;
end;
end

wyyang53
01-24-2008, 01:52 AM
As the ist author of 'Applied Numerical Methods Using MATLAB', I apologize for having left out the file 'grad.m'.
It contains the following statements:

function g= grad(f,x,h)
%grad.m to get the gradient of f(x) at x.
if nargin<3, h=.00001; end
h2= 2*h; N= length(x);
I= eye(N);
for n=1:N
g(n)= (feval(f,x+I(n,:)*h)-feval(f,x-I(n,:)*h))/h2;
end

Thanks.

cloyst3r
08-01-2008, 05:24 AM
What is behind the smilies in line g(n)= (feval(f,x+I(n............. ???
Can you rewrite this line WITHOUT smilies ???... Cuz it makes no sense replacing the smilies with : )

davekw7x
08-01-2008, 08:07 PM
What is behind the smilies...

For people trying to read code that was posted without code tags:
Pretend you are going to respond to the post (click the Quote button at the lower right).

When the response edit window opens, copy/paste the problematic lines to your text editor. Then, of course, unless you are really going to respond to the post, you click the Cancel button.

Here's what was in that post:




function g= grad(f,x,h)
%grad.m to get the gradient of f(x) at x.
if nargin<3, h=.00001; end
h2= 2*h; N= length(x);
I= eye(N);
for n=1:N
g(n)= (feval(f,x+I(n,:)*h)-feval(f,x-I(n,:)*h))/h2;
end


For people posting code: The board will preserve your indentation and the post will show up in a monospace font (like the stuff in the Code: box, above) if you highlight your lines of code (just the code; not your other narrative), and click the # icon at the top of the edit window. This puts code tags around the highlighted stuff.

When posting, If you click the Preview Post button before submitting the reply, you can see if you got it right.

Cuz it makes no sense replacing the smilies with : )Actually, it does make sense to Matlab, I'm thinking. I mean, in general, m(n,:) designates the nth row of matrix m, doesn't it?

Regards,

Dave

cloyst3r
08-02-2008, 01:39 AM
Thanks Dave ! ;)