Commit 09f1e806 authored by Ulmer Louis's avatar Ulmer Louis
Browse files

add OPTIM

parent 2373ff12
No preview for this file type
File added
{
"cells": [],
"metadata": {},
"nbformat": 4,
"nbformat_minor": 2
}
This source diff could not be displayed because it is too large. You can view the blob instead.
This diff is collapsed.
This diff is collapsed.
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# TP2 : Minimisation de la fonction log barrière"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Exercice 1 : Minimisation de la fonction log barrière"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"import cvxpy as cvx\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 1."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"ename": "NameError",
"evalue": "name 'p' is not defined",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m<ipython-input-3-be2d6f81c260>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0mc\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;36m6.556610413869322\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 2\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 3\u001b[0;31m \u001b[0malpha\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcvx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mVariable\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mp\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 4\u001b[0m \u001b[0mobj\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcvx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mMinimize\u001b[0m\u001b[0;34m(\u001b[0m \u001b[0;34m-\u001b[0m\u001b[0mc\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcvx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlog\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0mpow\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0malpha\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m/\u001b[0m\u001b[0mpow\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mc\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 5\u001b[0m \u001b[0mprob\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcvx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mProblem\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mNameError\u001b[0m: name 'p' is not defined"
]
}
],
"source": [
"c = 6.556610413869322\n",
"\n",
"alpha = cvx.Variable(p+1)\n",
"obj = cvx.Minimize( -c*sum(cvx.log(1-pow(X*alpha-y,2)/pow(c,2))) )\n",
"prob = cvx.Problem(obj)\n",
"prob.solve()\n",
"print(alpha.value)\n",
"\n",
"plt.plot(x,f)\n",
"plt.plot(x,y,'ro')\n",
"z = Xt*alpha.value\n",
"plt.plot(xt,z,'g')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 2."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def cout_logB(a,X,y,c):\n",
" r = X@a-y\n",
" cout = -c*np.sum(np.log(1-pow(r,2)/pow(c,2)))\n",
" return cout"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 3."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### a."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def grad_logB(a,X,y,c):\n",
" r = X@a-y;\n",
" df = 2*c*r/(pow(c,2)-pow(r,2));\n",
" g = X.T@df;\n",
" return g"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### b."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 4."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"def Hess_logB(a,X,y,c):\n",
" r = X@a-y;\n",
" d2f = 2*c*(pow(c,2)+pow(r,2))/(pow(pow(c,2)-pow(r,2),2));\n",
" D = np.diag(d2f);\n",
" H = X.T@D@X;\n",
" return H"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 5."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### a."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### b."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### c."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### d."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 6."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### a."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### b."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 7."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 8."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### 9."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"ename": "NameError",
"evalue": "name 'X' is not defined",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mNameError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m<ipython-input-3-7e4924456b9f>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m 18\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 19\u001b[0m \u001b[0malpha\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcvx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mVariable\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mp\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 20\u001b[0;31m \u001b[0mobj\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcvx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mMinimize\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mcvx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msum_squares\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mX\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0malpha\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 21\u001b[0m \u001b[0mprob\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcvx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mProblem\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 22\u001b[0m \u001b[0mprob\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msolve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mNameError\u001b[0m: name 'X' is not defined"
]
}
],
"source": [
"import cvxpy as cvx\n",
"import numpy as np\n",
"\n",
"np.random.seed(42)\n",
"\n",
"n = 50\n",
"p = 12\n",
"x = np.asarray(sorted(np.random.rand(n)))\n",
"e = np.random.randn(n)\n",
"f = np.cos(np.pi*x)\n",
"y = f + np.sqrt(0.01)*e\n",
"\n",
"X = np.vander(x, p+1)\n",
"\n",
"nt = 1000\n",
"xt = np.linspace(0,1,nt)\n",
"Xt = np.vander(xt, p+1)\n",
"\n",
"alpha = cvx.Variable(p+1)\n",
"obj = cvx.Minimize(cvx.sum_squares(X+alpha - y))\n",
"prob = cvx.Problem(obj)\n",
"prob.solve()\n",
"print(alpha.value)\n",
"\n",
"import matplotlib.pyplot as plt\n",
"\n",
"plt.plot(x,f)\n",
"plt.plot(x,y,'ro')\n",
"z = Xt*alpha.value\n",
"plt.plot(xt,z,'g')\n",
"\n",
"def cout_logB(a,X,y,c):\n",
" r = X@a-y\n",
" cout = -c*np.sum(np.log(1-pow(r,2)/pow(c,2)))\n",
" return cout\n",
"\n",
"def grad_logB(a,X,y,c):\n",
" r = X@a-y;\n",
" df = 2*c*r/(pow(c,2)-pow(r,2));\n",
" g = X.T@df;\n",
" return g\n",
"\n",
"def Hess_logB(a,X,y,c):\n",
" r = X@a-y;\n",
" d2f = 2*c*(pow(c,2)+pow(r,2))/(pow(pow(c,2)-pow(r,2),2));\n",
" D = np.diag(d2f);\n",
" H = X.T@D@X;\n",
" return H\n",
"\n",
"\n",
"a = np.asarray(alpha.value)\n",
"a = a.reshape((13,))\n",
"a = 2*a\n",
"\n",
"c = np.max(np.abs(X@a-y))+0.1\n",
"\n",
"#test grablogB\n",
"gr = grad_logB(a,X,y,c)\n",
"c1 = cout_logB(a,X,y,c)\n",
"eps= 0.0000001\n",
"e = 0*a\n",
"e[0]=1\n",
"c2 = cout_logB(a+eps*e,X,y,c)\n",
"\n",
"gr[0] - (c2-c1)/eps\n",
"\n",
"#\n",
"a = np.random.rand(p+1)\n",
"c = np.max(np.abs(X@a-y))+0.1\n",
"ainit = a;\n",
"\n",
"g = 1\n",
"k = 1;\n",
"nbitemax = 1000;\n",
"v,_ = np.linalg.eig(X.T@X)\n",
"pas = 1/np.max(np.abs(v))\n",
"pas = 0.001\n",
"\n",
"while ((np.linalg.norm(g) > 0.01) & (k < nbitemas)):\n",
" cout = cout_logB(a,X,y,c)\n",
" print(\"iteration %4d cout %7.4f\" % (k, cout))\n",
" g = grad_logB(a,X,y,c)\n",
" a = a - pas*g\n",
" k = k+1\n",
"\n",
"plt.plot(x,f)\n",
"plt.plot(x,y,'ro')\n",
"z = Xt@a\n",
"plt.plot(xt,z,'g')\n",
"plt.show()\n",
"\n",
"\n",
"\n",
"a = ainit\n",
"g = 1\n",
"k = 1\n",
"\n",
"alpha = .15\n",
"beta = 2\n",
"cout = cout_logB(a,X,y,c)\n",
"cout_old = cout+1\n",
"\n",
"while ((np.linalg.norm(g) > 0.01) & (k < nbitemas)):\n",
" cout = cout_logB(a,X,y,c)\n",
" print(\"iteration %4d cout %7.4f\" % (k, cout))\n",
" \n",
" if (cout < cout_old):\n",
" pas = (1+alpha)*pas\n",
" cout_old = cout\n",
" else:\n",
" a = a + pas*g\n",
" pas = pas/beta\n",
" \n",
" g = grad_logB(a,X,y,c)\n",
" a = a - pas*g\n",
" k = k+1\n",
" \n",
" \n",
"plt.plot(x,f)\n",
"plt.plot(x,y,'ro')\n",
"z = Xt@a\n",
"plt.plot(xt,z,'g')\n",
"plt.show()\n",
"\n",
"a = ainit\n",
"I = np.eye(p+1)\n",
"lamb = 0.001\n",
"k = 1\n",
"g = 2\n",
"\n",
"while ((np.linalg.norm(g) > 0.01) & (k < nbitemas)):\n",
" cout = cout_logB(a,X,y,c)\n",
" print(\"iteration %4d cout %7.4f\" % (k, cout))\n",
" g = grad_logB(a,X,y,c)\n",
" H = Hess_logB(a,X,y,c)\n",
" dir = np.linalg.solve(H+lamb*I,g)\n",
" a = a - dir\n",
" k = k+1\n",
" \n",
"plt.plot(x,f)\n",
"plt.plot(x,y,'ro')\n",
"z = Xt@a\n",
"plt.plot(xt,z,'g')\n",
"plt.show"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.5"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
%% Cell type:markdown id: tags:
# TP2 : Minimisation de la fonction log barrière
%% Cell type:markdown id: tags:
## Exercice 1 : Minimisation de la fonction log barrière
%% Cell type:code id: tags:
``` python
import cvxpy as cvx
import numpy as np
import matplotlib.pyplot as plt
```
%% Cell type:markdown id: tags:
### 1.
%% Cell type:code id: tags:
``` python
c = 6.556610413869322
alpha = cvx.Variable(p+1)
obj = cvx.Minimize( -c*sum(cvx.log(1-pow(X*alpha-y,2)/pow(c,2))) )
prob = cvx.Problem(obj)
prob.solve()
print(alpha.value)
plt.plot(x,f)
plt.plot(x,y,'ro')
z = Xt*alpha.value
plt.plot(xt,z,'g')
```
%%%% Output: error
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
<ipython-input-3-be2d6f81c260> in <module>()
1 c = 6.556610413869322
2
----> 3 alpha = cvx.Variable(p+1)
4 obj = cvx.Minimize( -c*sum(cvx.log(1-pow(X*alpha-y,2)/pow(c,2))) )
5 prob = cvx.Problem(obj)
NameError: name 'p' is not defined
%% Cell type:markdown id: tags:
### 2.
%% Cell type:code id: tags:
``` python
def cout_logB(a,X,y,c):
r = X@a-y
cout = -c*np.sum(np.log(1-pow(r,2)/pow(c,2)))
return cout
```
%% Cell type:markdown id: tags:
### 3.
%% Cell type:markdown id: tags:
#### a.
%% Cell type:code id: tags:
``` python
def grad_logB(a,X,y,c):
r = X@a-y;
df = 2*c*r/(pow(c,2)-pow(r,2));
g = X.T@df;
return g
```
%% Cell type:markdown id: tags:
#### b.
%% Cell type:code id: tags:
``` python
```
%% Cell type:markdown id: tags:
### 4.
%% Cell type:code id: tags:
``` python
def Hess_logB(a,X,y,c):
r = X@a-y;
d2f = 2*c*(pow(c,2)+pow(r,2))/(pow(pow(c,2)-pow(r,2),2));
D = np.diag(d2f);
H = X.T@D@X;
return H
```
%% Cell type:markdown id: tags:
### 5.
%% Cell type:markdown id: tags:
#### a.
%% Cell type:code id: tags:
``` python
```
%% Cell type:markdown id: tags:
#### b.
%% Cell type:code id: tags:
``` python
```
%% Cell type:markdown id: tags:
#### c.
%% Cell type:code id: tags:
``` python
```
%% Cell type:markdown id: tags:
#### d.
%% Cell type:code id: tags:
``` python
```
%% Cell type:markdown id: tags:
### 6.
%% Cell type:markdown id: tags:
#### a.
%% Cell type:code id: tags:
``` python
```
%% Cell type:markdown id: tags:
#### b.
%% Cell type:code id: tags:
``` python
```
%% Cell type:markdown id: tags:
### 7.
%% Cell type:code id: tags:
``` python
```
%% Cell type:markdown id: tags:
### 8.
%% Cell type:code id: tags:
``` python
```
%% Cell type:markdown id: tags:
### 9.
%% Cell type:code id: tags:
``` python
import cvxpy as cvx
import numpy as np
np.random.seed(42)
n = 50
p = 12
x = np.asarray(sorted(np.random.rand(n)))
e = np.random.randn(n)
f = np.cos(np.pi*x)
y = f + np.sqrt(0.01)*e
X = np.vander(x, p+1)
nt = 1000