Fortran Optimizing Steepest Descent in Fortran: How to Find the Optimal Step Size?

[massalami]

Hello All,

I'm new to the forum and new to fortran as well.
I've been trying to create a method for implementing the steepest descent method for optimization.
I've been able to create the code for a fixed step size and it worked beaufitully, but I can not seem to be able to find the optimal step size.
My teacher asked us to use a cubic interpolation for it, and I was able to create a system using a stepsize interval and finding the best in the interval, but in the program, the loop to find this step size crashes.

If needed, I can post the code to what I've done, but it's kind of messy( I am VERY unorganized). lol
Any tips?

 

[gsal]

Post the code. Use start CODE and end /CODE tags each surrounded by square brackets ('[' and ']').

 

[massalami]

ok then.
Basically explaining what i used, was to set alfa1 to 0 and alfa2 to the upper limit. Then I used a cubic interpolation, using a function of alfa1 and a function of alfa 2, with the same coeficients and a function of their derivative in respect to alfa 1 and alfa2.

the function of alfa1 is equal to the function using x-alfa1*grad and the derivative of this function is equal to the derivative of the function using x-alfa*grad.
The same goes for alfa2.

The criteria to stopping is when the derivative of the function with respect to alfa1<0 and for alfa2 is >0
the alfa for the general function then becomes alfa=(alfa1+alfa2)/2
and then it is used in the program(that part is ok).

Just help me out with the loop. It will make more sense when reading the code

program steepest_descent
IMPLICIT NONE

  integer k,kmax
  real(8) alfa1,alfa2,f,x(2),x0(2), grad(2), delta,v1,v2,v3,v4,a,b,c,d,alfa
  kmax=20
  k=0
  delta=0.00001
  alfa1=0
  write(*,*)'Entre x, y, alfa'
  read(*,*)x(1),x(2),alfa2
  do while(k<kmax)
    k=k+1
    x0=x
    call alpha1(x0,alfa1,v1,v2)
    call alpha2(x0,alfa2,v3,v4)
	do while (alfa1.ne.alfa2)		
			if (v4>0)then 
				if(v2<0)then
					alfa1=(alfa1+alfa2)/2
					alfa2=alfa1
				else
					alfa1=alfa1-0.01
				end if
			else
				alfa1=alfa2
				alfa2=alfa2+0.01
			end if
	 end do
	alfa=alfa1
	x=x0-(alfa*grad)
    call funcao(x,f)
	write(*,4)k,x(1),x(2),f
4   format(i2,1x,3(f15.8))
  end do
end program
subroutine gradx(x,grad)
  real(8) x(2),grad(2)
  grad=x**3-9.*x**2+22.*x-13.
end subroutine gradx
subroutine alpha1(alfa1,v1,v2)
	real x(2),v1,v2,alfa1,grad(2),x0(2)
	call gradx(x0,grad)
	x=x0-alfa1*grad
	call funcao(x,f)
	v1=f
	x=x0-(alfa1+0.00001)*grad
	call funcao(x,f)
	v2=(f-v1)/0.00001
end subroutine alpha1
subroutine alpha2(x0,alfa2,v3,v4)
	real x(2),v4,v3,alfa2,grad(2),x0(2)
	call gradx(x0,grad)
	x=x0-alfa2*grad
	call funcao(x,f)
	v3=f
	x=x0-(alfa2+0.00001)*grad
	call funcao(x,f)
	v4=(f-v3)/0.00001
end subroutine alpha2

subroutine variaveis(alfa1,alfa2,v1,v2,v3,v4,a,b,c,d)
  real(8) alfa1,alfa2,v1,v2,v3,v4,a,b,c,d,x,y
x=alfa1
y=alfa2
a1=-(x*(-v2*(y**3)-3.*v1*(y**2))+v1*(y**3)+(x**2)*(-v4*(y**2)+v2*(y**2)+3.*v3*y)+(x**3)*(v4*y-v3))
a2=(-(y**3)+3.*x*(y**2)-3.*(x**2)*y+(x**3))
a=a1/a2
b=(-v2*(y**3)+x*(-2.*v4*(y**2)-v2*(y**2)+6.*v3*y-6.*v1*y)+(x**2)*(v4*(y+2.*v2*y)+v4*(x**3)))/(-(y**3)+3.*x*(y**2)-3.*(x**2)*y+(x**3))
c=-(-v4*(y**2)-2.*v2*(y**2)+x*(-v4*y+v2*y+3.*v3-3.*v1)+3.*v3*y-3.*v1*y+(2.*v4+v2)*(x**2))/(-(y**3)+3.*x*(y**2)-3.*(x**2)*y+(x**3))
d=(-v4*y-v2*y+(v4+v2)*x+2.*v3-2.*v1)/(-(y**3)+3.*x*(y**2)-3.*(x**2)*y+(x**3))
end subroutine variaveis
subroutine funcao(x,f)
  real(8) x(2),f
  f=0.25*x(1)**4-3.*x(1)**3+11.*x(1)**2-13.*x(1)+  &
    0.25*x(2)**4-3.*x(2)**3+11.*x(2)**2-13.*x(2)
end subroutine funcao

 

[gsal]

are you doing any thing special to be able to go beyond 132 columns? Just checking, I did not think that happened by itself

 

[gsal]

the declaration of subroutine alpha1() is missing the first argument, you only have 3 but call it with 4

 

[gsal]

you need to be consistent in your declarations...either use real(8) everywhere or nowhere; otherwise you are mixing here and there

 

[gsal]

your variable grad in the main program is never getting set to anything

whatever modifications are happening to grad, they are staying local to alpha1() and alpha2() since you never passed grad to these subs...so, the variable grad in there is just a local one and does not refer to the one in the main program

 

[massalami]

I appreciate the help and comments.
I apologize for the mess. This program was kind of rushed and I wasn't able to develop it properly.

I will try to do as you guys told me and re do it, but otherwise, I will come and bother you guys again

 

[massalami]

Ok.
GSal, I did as you instructed me, corrected the formatting, fixed the alfa1 call to add 4 parameters but I am still getting an error.
The program stops when i run, but compiles fine, although I understand that is irrelevant.

I believe the error is in the do while loop that includes the alfa1 and alfa2, since when choosing a set alfa, instead of doing the adaptative step size, the program ran and worked fine.
Can you check that specific section to see what mistake i did there, please?

 

[uart]

do while (alfa1.ne.alfa2)

Hi massalami, comparing reals for equality can be problematic and can sometimes give unexpected results. I haven't read through your code in detail, but is that only meant to detect equality after the assignment, "alfa1 = alfa2", or is it meant to detect equality after more general mathematical operations. The latter may not work. Generally better to use something like "abs(alfa2-alfa1) > toln" rather than ".ne."

 

[massalami]

alright, uart.
thanks for the input.
I'll try fixing that. Now that I've retried the program, I'm still having trouble with the program without the loop, meaning that my alpha1 and alpha2 functions and variable functions might have some problems with it too.
Even without the if and the equality part(which I've tried without to try to identify the problem).
any ideas?
thank you all A LOT!

 

[gsal]

post latest code
post command used for compilation