#OK to post homework
#Yukun Yao, Mar. 8, Assignment 13



###Problem 1### This is an open-ended question. In class it was very disappointing that G1 and G2 did not do well. By starting with very simple functions, like (x-1)^2+(y-1)^2 of two variables (whose minimum is obviously [1,1]), and then moving on to more variables and more complicated polynomials, see if either G1 or G2 gets the answer, and what choices of sz and tol make it work. Perhaps there is a better way to implement this than either G1 or G2. You can try and read section 7.3 in the Machine learning book in the secret EM19 directory.

#I tried G2((x-1)^2+(y-1)^2+1, [x, y], [0, 0], 0.1e-1, 1000, 0.1e-7) and got [.9999999963, .9999999963] which is quite good. I believe we should let tol be very small since we still have MaxI to bound the total number of iterations. In general, G2 works in more cases. Small tol and sz around 0.01 usually are good. Of course, MaxI should be reasonably large. 














### From C13.txt ###


#C13.txt: March 7, 2019 Dr. Z.'s Math640(RU)
#Implementing the gradient descent algorithm
Help:=proc(): print(`OptF(f,var), Grad(f,var) , GD1(f,var,pt,sz), GD2(f,var,pt,sz) `): 
print(`G1(f,var,pt,sz,MaxI,tol), G2(f,var,pt,sz,MaxI,tol) `):
end:

#OptF(f,var): inputs a nice function f (expression)
# of a list of variables, and finds the local max and min locations
#and values. For example (x^2+y^2+1,[x,y]); should output {[[0,0],1]}.
OptF:=proc(f,var) local eq,i,sol,L:
eq:={seq(diff(f,var[i])=0,i=1..nops(var))}:
sol:=solve(eq,var):

if sol=NULL then
  RETURN(FAIL):
fi:

L:=[seq(subs(sol[i],var),i=1..nops(sol))]: 
L:
end:

Grad:=proc(f,var) local i: option remember:
 [seq(diff(f,var[i]),i=1..nops(var))]: end:

#GD1(f,var,pt,sz): inputs a function (expression) f in the list of variables var, a point pt
#a step-size sz, and a small parameter e indicating when to stop
#Performs one step, returns the next point NORMALIZED
GD1:=proc(f,var,pt,sz) local g,i,Ng:
if not ( type(var,list) and type(pt,list) and nops(var)=nops(pt) and type(sz,numeric) and sz>0) then
  print(`bad input`):
  RETURN(FAIL):
fi:
  
 g:=Grad(f,var):

 g:=subs({seq(var[i]=pt[i],i=1..nops(var))},g):
 Ng:=evalf(sqrt(add(g[i]^2,i=1..nops(g)))):
 g:=g/Ng:

 pt-sz*g:

 end:

#GD2(f,var,pt,sz): inputs a function (expression) f in the list of variables var, a point pt
#a step-size sz, and a small parameter e indicating when to stop
#Performs one step, returns the next point NOT Normalized
GD2:=proc(f,var,pt,sz) local g,i:
if not ( type(var,list) and type(pt,list) and nops(var)=nops(pt) and type(sz,numeric) and sz>0) then
  print(`bad input`):
  RETURN(FAIL):
fi:
  
 g:=Grad(f,var):

 g:=subs({seq(var[i]=pt[i],i=1..nops(var))},g):
 

 pt-sz*g:

 end:

    #G1(f,var,pt,sz,MaxI,tol): inputs a function f given as an expression, in a list of variables var
    #a starting point pt, MaxI: number of iterations, and tol a small parameter, we stop when the
    #norm of the gradient is less than tol, or hit MaxI, iterations. NORMALIZED
   G1:=proc(f,var,pt,sz,MaxI,tol) local g, g1, pt1,i,N1:
   pt1:=pt:
     g:=Grad(f,var):
       g1:=subs({seq(var[i]=pt1[i],i=1..nops(var))},g):
       g1:=evalf(sqrt(add(g1[i]^2,i=1..nops(g1)))):
     for i from 1 to MaxI while g1>tol do
       pt1:=GD1(f,var,pt1,sz):
       
       g1:=subs({seq(var[i]=pt1[i],i=1..nops(var))},g):
       g1:=evalf(sqrt(add(g1[i]^2,i=1..nops(g1)))):
   od:
   pt1:    
          end:
 #G2(f,var,pt,sz,MaxI,tol): inputs a function f given as an expression, in a list of variables var
    #a starting point pt, MaxI: number of iterations, and tol a small parameter, we stop when the
    #norm of the gradient is less than tol, or hit MaxI, iterations. NOT NORMALIZED
   G2:=proc(f,var,pt,sz,MaxI,tol) local g, g1, pt1,i,N1:
   pt1:=pt:
     g:=Grad(f,var):
       g1:=subs({seq(var[i]=pt1[i],i=1..nops(var))},g):
       g1:=evalf(sqrt(add(g1[i]^2,i=1..nops(g1)))):
     for i from 1 to MaxI while g1>tol do
       pt1:=GD2(f,var,pt1,sz):
       
       g1:=subs({seq(var[i]=pt1[i],i=1..nops(var))},g):
       g1:=evalf(sqrt(add(g1[i]^2,i=1..nops(g1)))):
   od:
   pt1:    
          end: