Gradient descent: Difference between revisions

Line 534:

=={{header|Phix}}==

<lang Phix>-- Function for which minimum is to be found.

with javascript_semantics

function g(sequence x)

-- Function for which minimum is to be found.

atom {x0,x1} = x

function g(sequence x)

return (x0-1)*(x0-1)*exp(-x1*x1) +

atom {x0,x1} = x

x1*(x1+2)*exp(-2*x0*x0)

return (x0-1)*(x0-1)*exp(-x1*x1) +

end function

x1*(x1+2)*exp(-2*x0*x0)

end function

-- Provides a rough calculation of gradient g(x).

function gradG(sequence p)

-- Provides a rough calculation of gradient g(x).

atom {x,y} = p

function gradG(sequence p)

p[1] = 2*(x-1)*exp(-y*y) - 4*x*exp(-2*x*x)*y*(y+2)

atom {x,y} = p,

p[2] = -2*(x-1)*(x-1)*y*exp(-y*y) + exp(-2*x*x)*(y+2) + exp(-2*x*x)*y

xm1 = x-1,

return p

emyy = exp(-y*y),

end function

em2xx = exp(-2*x*x),

xm12emyy = xm1*2*emyy

function steepestDescent(sequence x, atom alpha, tolerance)

p = {xm12emyy - 4*x*em2xx*y*(y+2),

integer n = length(x)

-xm12emyy*xm1*y + em2xx*(2*y+2)}

atom g0 = g(x) -- Initial estimate of result.

return p

end function

-- Calculate initial gradient.

sequence fi = gradG(x)

function steepestDescent(sequence x, atom alpha, tolerance)

sequence fi = gradG(x) -- Calculate initial gradient

-- Calculate initial norm.

atom g0 = g(x), -- Initial estimate of result

atom delG = sqrt(sum(sq_mul(fi,fi))),

delG = sqrt(sum(sq_mul(fi,fi))), -- & norm

b = alpha / delG

b = alpha / delG

integer icount = 0 -- iteration limitor/sanity

-- Iterate until value is <= tolerance.

while delG>tolerance do -- Iterate until <= tolerance

while delG>tolerance do

x = sq_sub(x,sq_mul(b,fi)) -- Calculate next value

-- Calculate next value.

fi = gradG(x) -- next gradient

x = sq_sub(x,sq_mul(b,fi))

delG = sqrt(sum(sq_mul(fi,fi))) -- next norm

b = alpha / delG

-- Calculate next gradient.

atom g1 = g(x) -- next value

fi = gradG(x)

if g1>g0 then

alpha /= 2 -- Adjust parameter

-- Calculate next norm.

else

delG = sqrt(sum(sq_mul(fi,fi)))

g0 = g1

b = alpha / delG

end if

icount += 1

-- Calculate next value.

assert(icount<100) -- (increase if/when necessary)

atom g1 = g(x)

end while

return x

-- Adjust parameter.

end function

if g1>g0 then

alpha /= 2

constant alpha = 0.1, tolerance = 0.00000000000001

else

sequence x = steepestDescent({0.1,-1}, alpha, tolerance)

g0 = g1

printf(1,"Testing steepest descent method:\n")

end if

printf(1,"The minimum is at x = %.13f, y = %.13f for which f(x, y) = %.15f\n", {x[1], x[2], g(x)})

end while

return x

end function

constant tolerance = 0.0000001, alpha = 0.1

sequence x = steepestDescent({0.1,-1}, alpha, tolerance)

printf(1,"Testing steepest descent method:\n")

printf(1,"The minimum is at x = %.13f, y = %.13f for which f(x, y) = %.16f\n", {x[1], x[2], g(x)})</lang>

Results ~~now~~ match ~~(at least) Algol 68/W,~~ Fortran, ~~Go,~~ ~~Julia,~~ ~~Raku, REXX, and Wren [~~to ~~6dp~~ or ~~better anyway].~~

Results match Fortran, most others to 6 or 7dp

Some slightly unexpected/unusual (but I think acceptable) variance was noted when playing with different tolerances, be warned.

Note that specifying a tolerance < 1e-7 causes an infinite loop on Phix, whereas REXX copes with a much smaller tolerance.

Results on 32/64 bit Phix agree to 13dp, which I therefore choose to show in full here (but otherwise would not really trust).

<pre>

Testing steepest descent method:

The minimum is at x = 0.~~1076268243295~~, y = -1.~~2232596548816~~ for which f(x, y) = -0.~~7500634205514924~~

The minimum is at x = 0.1076268435484, y = -1.2232596638399 for which f(x, y) = -0.750063420551493

</pre>