foo

2cfefa0e · J. Gergaud · 0a9aa66a · 2cfefa0e
Commit 2cfefa0e authored 8 months ago by J. Gergaud
--- a/M2/TP-Newton.ipynb
+++ b/M2/TP-Newton.ipynb
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Automatic differentiation"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "using LinearAlgebra\n",
+    "using ForwardDiff\n",
+    "A = [-1 2;1 4]; b=[1,1]; c=1\n",
+    "f(x) = 0.5*x'*A*x + b'*x +c\n",
+    "\n",
+    "analytic_∇f(x) = 0.5*(A' + A)*x+b\n",
+    "analytic_∇²f(x) = 0.5*(A' + A)\n",
+    "\n",
+    "∇f(x) = ForwardDiff.gradient(f, x)\n",
+    "∇²f(x) = ForwardDiff.hessian(f, x)\n",
+    "\n",
+    "x0 = [1,-1]\n",
+    "println(\"∇f(x0) = \", ∇f(x0))\n",
+    "println(\"analytic_∇f(x0) = \", analytic_∇f(x0))\n",
+    "println(\"analytic_∇f(x0) - ∇f(x0) = \", analytic_∇f(x0)- ∇f(x0))\n",
+    "\n",
+    "println(\"∇²f(x0) = \", ∇²f(x0))\n",
+    "println(\"analytic_∇²f(x0) = \", analytic_∇²f(x0))\n",
+    "println(\"analytic_∇²f(x0) - ∇²f(x0) = \", analytic_∇²f(x0)- ∇²f(x0))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Newton algorithm"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "using LinearAlgebra\n",
+    "using ForwardDiff\n",
+    "\n",
+    "\"\"\"\n",
+    "   Solve by Newton's algorithm the optimization problem Min f(x)\n",
+    "   Case where f is a function from R^n to R\n",
+    "\"\"\"\n",
+    "\n",
+    "\n",
+    "function algo_Newton(f,x0::Vector{<:Real};AbsTol= abs(eps()), RelTol = abs(eps()), ε=0.01, nbit_max = 0)\n",
+    "# to complete\n",
+    "    \n",
+    "    # flag = 0 if the program stop on the first criteria\n",
+    "    # flag = 2 if the program stop on the second criteria\n",
+    "    # ...\n",
+    "    return xₖ, flag, fₖ, ∇fₖ , k  \n",
+    "end"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "A = [1 0 ; 0 9/2]\n",
+    "b = [0,0]; c=0.\n",
+    "f1(x) = 0.5*x'*A*x + b'*x +c\n",
+    "x0 = [1000,-20]\n",
+    "println(\"Results for Newton on f1 : \", algo_Newton(f1,x0))\n",
+    "println(\"eigen value of 0.5(A^T+A) = \", 0.5*eigen(A'+A).values)\n",
+    "using Plots;\n",
+    "x=range(-10,stop=10,length=100)\n",
+    "y=range(-10,stop=10,length=100)\n",
+    "f11(x,y) = f1([x,y])\n",
+    "p1 = plot(x,y,f11,st=:contourf)\n",
+    "\n",
+    "A = [-1 0 ; 0 3]\n",
+    "f2(x) = 0.5*x'*A*x + b'*x +c\n",
+    "println(\"Results for Newton on f2 : \", algo_Newton(f2,x0))\n",
+    "println(\"eigen value of 0.5(A^T+A) = \", 0.5*eigen(A'+A).values)\n",
+    "f21(x,y) = f1([x,y])\n",
+    "p2 = plot(x,y,f21,st=:contourf)\n",
+    "# Rosenbrock function\n",
+    "x=range(-1.5,stop=2,length=100)\n",
+    "y=range(-0.5,stop=3.,length=100)\n",
+    "f3(x) = (1-x[1])^2 + 100*(x[2]-x[1]^2)^2\n",
+    "f31(x,y) = f3([x,y])\n",
+    "p3 = plot(x,y,f31,st=:contourf)\n",
+    "x0 = [1.1,1.2]\n",
+    "\n",
+    "println(\"Results for Newton on f3 : \", algo_Newton(f3,[1.1,1.2]))\n",
+    "\n",
+    "println(\"Results for Newton on f3 : \", algo_Newton(f3,[3.,0.]))\n",
+    "plot(p1,p2,p3)"
+   ]
+  }
+ ],
+ "metadata": {
+  "language_info": {
+   "name": "python"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
+%% Cell type:markdown id: tags:
+
+## Automatic differentiation
+
+%% Cell type:code id: tags:
+
+``` 
+using LinearAlgebra
+using ForwardDiff
+A = [-1 2;1 4]; b=[1,1]; c=1
+f(x) = 0.5*x'*A*x + b'*x +c
+
+analytic_∇f(x) = 0.5*(A' + A)*x+b
+analytic_∇²f(x) = 0.5*(A' + A)
+
+∇f(x) = ForwardDiff.gradient(f, x)
+∇²f(x) = ForwardDiff.hessian(f, x)
+
+x0 = [1,-1]
+println("∇f(x0) = ", ∇f(x0))
+println("analytic_∇f(x0) = ", analytic_∇f(x0))
+println("analytic_∇f(x0) - ∇f(x0) = ", analytic_∇f(x0)- ∇f(x0))
+
+println("∇²f(x0) = ", ∇²f(x0))
+println("analytic_∇²f(x0) = ", analytic_∇²f(x0))
+println("analytic_∇²f(x0) - ∇²f(x0) = ", analytic_∇²f(x0)- ∇²f(x0))
+```
+
+%% Cell type:markdown id: tags:
+
+## Newton algorithm
+
+%% Cell type:code id: tags:
+
+``` 
+using LinearAlgebra
+using ForwardDiff
+
+"""
+   Solve by Newton's algorithm the optimization problem Min f(x)
+   Case where f is a function from R^n to R
+"""
+
+
+function algo_Newton(f,x0::Vector{<:Real};AbsTol= abs(eps()), RelTol = abs(eps()), ε=0.01, nbit_max = 0)
+# to complete
+
+    # flag = 0 if the program stop on the first criteria
+    # flag = 2 if the program stop on the second criteria
+    # ...
+    return xₖ, flag, fₖ, ∇fₖ , k
+end
+```
+
+%% Cell type:code id: tags:
+
+``` 
+A = [1 0 ; 0 9/2]
+b = [0,0]; c=0.
+f1(x) = 0.5*x'*A*x + b'*x +c
+x0 = [1000,-20]
+println("Results for Newton on f1 : ", algo_Newton(f1,x0))
+println("eigen value of 0.5(A^T+A) = ", 0.5*eigen(A'+A).values)
+using Plots;
+x=range(-10,stop=10,length=100)
+y=range(-10,stop=10,length=100)
+f11(x,y) = f1([x,y])
+p1 = plot(x,y,f11,st=:contourf)
+
+A = [-1 0 ; 0 3]
+f2(x) = 0.5*x'*A*x + b'*x +c
+println("Results for Newton on f2 : ", algo_Newton(f2,x0))
+println("eigen value of 0.5(A^T+A) = ", 0.5*eigen(A'+A).values)
+f21(x,y) = f1([x,y])
+p2 = plot(x,y,f21,st=:contourf)
+# Rosenbrock function
+x=range(-1.5,stop=2,length=100)
+y=range(-0.5,stop=3.,length=100)
+f3(x) = (1-x[1])^2 + 100*(x[2]-x[1]^2)^2
+f31(x,y) = f3([x,y])
+p3 = plot(x,y,f31,st=:contourf)
+x0 = [1.1,1.2]
+
+println("Results for Newton on f3 : ", algo_Newton(f3,[1.1,1.2]))
+
+println("Results for Newton on f3 : ", algo_Newton(f3,[3.,0.]))
+plot(p1,p2,p3)
+```