small update on grad descent

This commit is contained in:
mhjensen
2019-09-20 06:39:17 +02:00
parent f6049d88b2
commit 9ba0c3c0d7
13 changed files with 174 additions and 431 deletions
File diff suppressed because one or more lines are too long
+1 -1
View File
@@ -292,7 +292,7 @@ MathJax.Hub.Config({
<center>[2] <b>Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University</b></center>
<br>
<p>
<center><h4>Sep 19, 2019</h4></center> <!-- date -->
<center><h4>Sep 20, 2019</h4></center> <!-- date -->
<br>
<p>
-21
View File
@@ -289,28 +289,7 @@ when \( ||\nabla_\beta C(\beta_k) || \leq \epsilon = 10^{-8} \).
<p>
And finally we can compare our solution for \( \beta \) with the analytic result given by
\( \beta= (X^TX)^{-1} X^T \mathbf{y} \).
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span><span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">numpy</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">np</span>
<span style="color: #BA2121; font-style: italic">&quot;&quot;&quot;</span>
<span style="color: #BA2121; font-style: italic">The following setup is just a suggestion, feel free to write it the way you like.</span>
<span style="color: #BA2121; font-style: italic">&quot;&quot;&quot;</span>
<span style="color: #408080; font-style: italic">#Setup problem described in the exercise</span>
N <span style="color: #666666">=</span> <span style="color: #666666">100</span> <span style="color: #408080; font-style: italic">#Nr of datapoints</span>
M <span style="color: #666666">=</span> <span style="color: #666666">2</span> <span style="color: #408080; font-style: italic">#Nr of features</span>
x <span style="color: #666666">=</span> np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>rand(N) <span style="color: #408080; font-style: italic">#Uniformly generated x-values in [0,1]</span>
y <span style="color: #666666">=</span> <span style="color: #666666">5*</span>x<span style="color: #666666">**2</span> <span style="color: #666666">+</span> <span style="color: #666666">0.1*</span>np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>randn(N)
X <span style="color: #666666">=</span> np<span style="color: #666666">.</span>c_[np<span style="color: #666666">.</span>ones(N),x] <span style="color: #408080; font-style: italic">#Construct design matrix</span>
<span style="color: #408080; font-style: italic">#Compute beta according to normal equations to compare with GD solution</span>
Xt_X_inv <span style="color: #666666">=</span> np<span style="color: #666666">.</span>linalg<span style="color: #666666">.</span>inv(np<span style="color: #666666">.</span>dot(X<span style="color: #666666">.</span>T,X))
Xt_y <span style="color: #666666">=</span> np<span style="color: #666666">.</span>dot(X<span style="color: #666666">.</span>transpose(),y)
beta_NE <span style="color: #666666">=</span> np<span style="color: #666666">.</span>dot(Xt_X_inv,Xt_y)
<span style="color: #008000; font-weight: bold">print</span>(beta_NE)
</pre></div>
<p>
<p>
<!-- navigation buttons at the bottom of the page -->
+5 -4
View File
@@ -289,17 +289,18 @@ Another simple example is here
<span style="color: #008000; font-weight: bold">from</span> <span style="color: #0000FF; font-weight: bold">matplotlib.ticker</span> <span style="color: #008000; font-weight: bold">import</span> LinearLocator, FormatStrFormatter
<span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">sys</span>
x <span style="color: #666666">=</span> <span style="color: #666666">2*</span>np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>rand(<span style="color: #666666">100</span>,<span style="color: #666666">1</span>)
y <span style="color: #666666">=</span> <span style="color: #666666">4+3*</span>x<span style="color: #666666">+</span>np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>randn(<span style="color: #666666">100</span>,<span style="color: #666666">1</span>)
<span style="color: #408080; font-style: italic"># the number of datapoints</span>
m <span style="color: #666666">=</span> <span style="color: #666666">100</span>
x <span style="color: #666666">=</span> <span style="color: #666666">2*</span>np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>rand(m,<span style="color: #666666">1</span>)
y <span style="color: #666666">=</span> <span style="color: #666666">4+3*</span>x<span style="color: #666666">+</span>np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>randn(m,<span style="color: #666666">1</span>)
xb <span style="color: #666666">=</span> np<span style="color: #666666">.</span>c_[np<span style="color: #666666">.</span>ones((<span style="color: #666666">100</span>,<span style="color: #666666">1</span>)), x]
xb <span style="color: #666666">=</span> np<span style="color: #666666">.</span>c_[np<span style="color: #666666">.</span>ones((m,<span style="color: #666666">1</span>)), x]
beta_linreg <span style="color: #666666">=</span> np<span style="color: #666666">.</span>linalg<span style="color: #666666">.</span>inv(xb<span style="color: #666666">.</span>T<span style="color: #666666">.</span>dot(xb))<span style="color: #666666">.</span>dot(xb<span style="color: #666666">.</span>T)<span style="color: #666666">.</span>dot(y)
<span style="color: #008000; font-weight: bold">print</span>(beta_linreg)
beta <span style="color: #666666">=</span> np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>randn(<span style="color: #666666">2</span>,<span style="color: #666666">1</span>)
eta <span style="color: #666666">=</span> <span style="color: #666666">0.1</span>
Niterations <span style="color: #666666">=</span> <span style="color: #666666">1000</span>
m <span style="color: #666666">=</span> <span style="color: #666666">100</span>
<span style="color: #008000; font-weight: bold">for</span> <span style="color: #008000">iter</span> <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #008000">range</span>(Niterations):
gradients <span style="color: #666666">=</span> <span style="color: #666666">2.0/</span>m<span style="color: #666666">*</span>xb<span style="color: #666666">.</span>T<span style="color: #666666">.</span>dot(xb<span style="color: #666666">.</span>dot(beta)<span style="color: #666666">-</span>y)
+2 -34
View File
@@ -290,43 +290,11 @@ $$
$$
<p>
We can now extend our program to minimize \( C_{\text{ridge}}(\beta) \) using gradient descent and compare with the analytical solution given by
We can easily extend our program to minimize \( C_{\text{ridge}}(\beta) \) using gradient descent and compare with the analytical solution given by
$$
\beta_{\text{ridge}} = \left(X^T X + \lambda I_{2 \times 2} \right)^{-1} X^T \mathbf{y},
\beta_{\text{ridge}} = \left(X^T X + \lambda I_{2 \times 2} \right)^{-1} X^T \mathbf{y}.
$$
for \( \lambda = {0,1,10,50,100} \) (\( \lambda = 0 \) corresponds to ordinary least squares).
We can then compute \( ||\beta_{\text{ridge}}|| \) for each \( \lambda \).
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span><span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">numpy</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">np</span>
<span style="color: #BA2121; font-style: italic">&quot;&quot;&quot;</span>
<span style="color: #BA2121; font-style: italic">The following setup is just a suggestion, feel free to write it the way you like.</span>
<span style="color: #BA2121; font-style: italic">&quot;&quot;&quot;</span>
<span style="color: #408080; font-style: italic">#Setup problem described in the exercise</span>
N <span style="color: #666666">=</span> <span style="color: #666666">100</span> <span style="color: #408080; font-style: italic">#Nr of datapoints</span>
M <span style="color: #666666">=</span> <span style="color: #666666">2</span> <span style="color: #408080; font-style: italic">#Nr of features</span>
x <span style="color: #666666">=</span> np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>rand(N)
y <span style="color: #666666">=</span> <span style="color: #666666">5*</span>x<span style="color: #666666">**2</span> <span style="color: #666666">+</span> <span style="color: #666666">0.1*</span>np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>randn(N)
<span style="color: #408080; font-style: italic">#Compute analytic beta for Ridge regression </span>
X <span style="color: #666666">=</span> np<span style="color: #666666">.</span>c_[np<span style="color: #666666">.</span>ones(N),x]
XT_X <span style="color: #666666">=</span> np<span style="color: #666666">.</span>dot(X<span style="color: #666666">.</span>T,X)
l <span style="color: #666666">=</span> <span style="color: #666666">0.1</span> <span style="color: #408080; font-style: italic">#Ridge parameter lambda</span>
Id <span style="color: #666666">=</span> np<span style="color: #666666">.</span>eye(XT_X<span style="color: #666666">.</span>shape[<span style="color: #666666">0</span>])
Z <span style="color: #666666">=</span> np<span style="color: #666666">.</span>linalg<span style="color: #666666">.</span>inv(XT_X<span style="color: #666666">+</span>l<span style="color: #666666">*</span>Id)
beta_ridge <span style="color: #666666">=</span> np<span style="color: #666666">.</span>dot(Z,np<span style="color: #666666">.</span>dot(X<span style="color: #666666">.</span>T,y))
<span style="color: #008000; font-weight: bold">print</span>(beta_ridge)
<span style="color: #008000; font-weight: bold">print</span>(np<span style="color: #666666">.</span>linalg<span style="color: #666666">.</span>norm(beta_ridge)) <span style="color: #408080; font-style: italic">#||beta||</span>
</pre></div>
<p>
<p>
<!-- navigation buttons at the bottom of the page -->
+1 -1
View File
@@ -292,7 +292,7 @@ MathJax.Hub.Config({
<center>[2] <b>Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University</b></center>
<br>
<p>
<center><h4>Sep 19, 2019</h4></center> <!-- date -->
<center><h4>Sep 20, 2019</h4></center> <!-- date -->
<br>
<p>
+8 -62
View File
@@ -148,7 +148,7 @@ MathJax.Hub.Config({
<center>[2] <b>Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University</b></center>
<br>
<p>&nbsp;<br>
<center><h4>Sep 19, 2019</h4></center> <!-- date -->
<center><h4>Sep 20, 2019</h4></center> <!-- date -->
<br>
<p>
@@ -747,28 +747,6 @@ when \( ||\nabla_\beta C(\beta_k) || \leq \epsilon = 10^{-8} \).
<p>
And finally we can compare our solution for \( \beta \) with the analytic result given by
\( \beta= (X^TX)^{-1} X^T \mathbf{y} \).
<p>
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
<div class="highlight" style="background: #eeeedd"><pre style="font-size: 80%; line-height: 125%"><span></span><span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
<span style="color: #CD5555">&quot;&quot;&quot;</span>
<span style="color: #CD5555">The following setup is just a suggestion, feel free to write it the way you like.</span>
<span style="color: #CD5555">&quot;&quot;&quot;</span>
<span style="color: #228B22">#Setup problem described in the exercise</span>
N = <span style="color: #B452CD">100</span> <span style="color: #228B22">#Nr of datapoints</span>
M = <span style="color: #B452CD">2</span> <span style="color: #228B22">#Nr of features</span>
x = np.random.rand(N) <span style="color: #228B22">#Uniformly generated x-values in [0,1]</span>
y = <span style="color: #B452CD">5</span>*x**<span style="color: #B452CD">2</span> + <span style="color: #B452CD">0.1</span>*np.random.randn(N)
X = np.c_[np.ones(N),x] <span style="color: #228B22">#Construct design matrix</span>
<span style="color: #228B22">#Compute beta according to normal equations to compare with GD solution</span>
Xt_X_inv = np.linalg.inv(np.dot(X.T,X))
Xt_y = np.dot(X.transpose(),y)
beta_NE = np.dot(Xt_X_inv,Xt_y)
<span style="color: #8B008B; font-weight: bold">print</span>(beta_NE)
</pre></div>
</section>
@@ -789,17 +767,18 @@ Another simple example is here
<span style="color: #8B008B; font-weight: bold">from</span> <span style="color: #008b45; text-decoration: underline">matplotlib.ticker</span> <span style="color: #8B008B; font-weight: bold">import</span> LinearLocator, FormatStrFormatter
<span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">sys</span>
x = <span style="color: #B452CD">2</span>*np.random.rand(<span style="color: #B452CD">100</span>,<span style="color: #B452CD">1</span>)
y = <span style="color: #B452CD">4</span>+<span style="color: #B452CD">3</span>*x+np.random.randn(<span style="color: #B452CD">100</span>,<span style="color: #B452CD">1</span>)
<span style="color: #228B22"># the number of datapoints</span>
m = <span style="color: #B452CD">100</span>
x = <span style="color: #B452CD">2</span>*np.random.rand(m,<span style="color: #B452CD">1</span>)
y = <span style="color: #B452CD">4</span>+<span style="color: #B452CD">3</span>*x+np.random.randn(m,<span style="color: #B452CD">1</span>)
xb = np.c_[np.ones((<span style="color: #B452CD">100</span>,<span style="color: #B452CD">1</span>)), x]
xb = np.c_[np.ones((m,<span style="color: #B452CD">1</span>)), x]
beta_linreg = np.linalg.inv(xb.T.dot(xb)).dot(xb.T).dot(y)
<span style="color: #8B008B; font-weight: bold">print</span>(beta_linreg)
beta = np.random.randn(<span style="color: #B452CD">2</span>,<span style="color: #B452CD">1</span>)
eta = <span style="color: #B452CD">0.1</span>
Niterations = <span style="color: #B452CD">1000</span>
m = <span style="color: #B452CD">100</span>
<span style="color: #8B008B; font-weight: bold">for</span> <span style="color: #658b00">iter</span> <span style="color: #8B008B">in</span> <span style="color: #658b00">range</span>(Niterations):
gradients = <span style="color: #B452CD">2.0</span>/m*xb.T.dot(xb.dot(beta)-y)
@@ -869,45 +848,12 @@ $$
<p>&nbsp;<br>
<p>
We can now extend our program to minimize \( C_{\text{ridge}}(\beta) \) using gradient descent and compare with the analytical solution given by
We can easily extend our program to minimize \( C_{\text{ridge}}(\beta) \) using gradient descent and compare with the analytical solution given by
<p>&nbsp;<br>
$$
\beta_{\text{ridge}} = \left(X^T X + \lambda I_{2 \times 2} \right)^{-1} X^T \mathbf{y},
\beta_{\text{ridge}} = \left(X^T X + \lambda I_{2 \times 2} \right)^{-1} X^T \mathbf{y}.
$$
<p>&nbsp;<br>
for \( \lambda = {0,1,10,50,100} \) (\( \lambda = 0 \) corresponds to ordinary least squares).
We can then compute \( ||\beta_{\text{ridge}}|| \) for each \( \lambda \).
<p>
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
<div class="highlight" style="background: #eeeedd"><pre style="font-size: 80%; line-height: 125%"><span></span><span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
<span style="color: #CD5555">&quot;&quot;&quot;</span>
<span style="color: #CD5555">The following setup is just a suggestion, feel free to write it the way you like.</span>
<span style="color: #CD5555">&quot;&quot;&quot;</span>
<span style="color: #228B22">#Setup problem described in the exercise</span>
N = <span style="color: #B452CD">100</span> <span style="color: #228B22">#Nr of datapoints</span>
M = <span style="color: #B452CD">2</span> <span style="color: #228B22">#Nr of features</span>
x = np.random.rand(N)
y = <span style="color: #B452CD">5</span>*x**<span style="color: #B452CD">2</span> + <span style="color: #B452CD">0.1</span>*np.random.randn(N)
<span style="color: #228B22">#Compute analytic beta for Ridge regression </span>
X = np.c_[np.ones(N),x]
XT_X = np.dot(X.T,X)
l = <span style="color: #B452CD">0.1</span> <span style="color: #228B22">#Ridge parameter lambda</span>
Id = np.eye(XT_X.shape[<span style="color: #B452CD">0</span>])
Z = np.linalg.inv(XT_X+l*Id)
beta_ridge = np.dot(Z,np.dot(X.T,y))
<span style="color: #8B008B; font-weight: bold">print</span>(beta_ridge)
<span style="color: #8B008B; font-weight: bold">print</span>(np.linalg.norm(beta_ridge)) <span style="color: #228B22">#||beta||</span>
</pre></div>
</section>
+8 -60
View File
@@ -210,7 +210,7 @@ MathJax.Hub.Config({
<center>[2] <b>Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University</b></center>
<br>
<p>
<center><h4>Sep 19, 2019</h4></center> <!-- date -->
<center><h4>Sep 20, 2019</h4></center> <!-- date -->
<br>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
@@ -753,28 +753,7 @@ when \( ||\nabla_\beta C(\beta_k) || \leq \epsilon = 10^{-8} \).
<p>
And finally we can compare our solution for \( \beta \) with the analytic result given by
\( \beta= (X^TX)^{-1} X^T \mathbf{y} \).
<p>
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
<div class="highlight" style="background: #eeeedd"><pre style="line-height: 125%"><span></span><span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
<span style="color: #CD5555">&quot;&quot;&quot;</span>
<span style="color: #CD5555">The following setup is just a suggestion, feel free to write it the way you like.</span>
<span style="color: #CD5555">&quot;&quot;&quot;</span>
<span style="color: #228B22">#Setup problem described in the exercise</span>
N = <span style="color: #B452CD">100</span> <span style="color: #228B22">#Nr of datapoints</span>
M = <span style="color: #B452CD">2</span> <span style="color: #228B22">#Nr of features</span>
x = np.random.rand(N) <span style="color: #228B22">#Uniformly generated x-values in [0,1]</span>
y = <span style="color: #B452CD">5</span>*x**<span style="color: #B452CD">2</span> + <span style="color: #B452CD">0.1</span>*np.random.randn(N)
X = np.c_[np.ones(N),x] <span style="color: #228B22">#Construct design matrix</span>
<span style="color: #228B22">#Compute beta according to normal equations to compare with GD solution</span>
Xt_X_inv = np.linalg.inv(np.dot(X.T,X))
Xt_y = np.dot(X.transpose(),y)
beta_NE = np.dot(Xt_X_inv,Xt_y)
<span style="color: #8B008B; font-weight: bold">print</span>(beta_NE)
</pre></div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
@@ -794,17 +773,18 @@ Another simple example is here
<span style="color: #8B008B; font-weight: bold">from</span> <span style="color: #008b45; text-decoration: underline">matplotlib.ticker</span> <span style="color: #8B008B; font-weight: bold">import</span> LinearLocator, FormatStrFormatter
<span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">sys</span>
x = <span style="color: #B452CD">2</span>*np.random.rand(<span style="color: #B452CD">100</span>,<span style="color: #B452CD">1</span>)
y = <span style="color: #B452CD">4</span>+<span style="color: #B452CD">3</span>*x+np.random.randn(<span style="color: #B452CD">100</span>,<span style="color: #B452CD">1</span>)
<span style="color: #228B22"># the number of datapoints</span>
m = <span style="color: #B452CD">100</span>
x = <span style="color: #B452CD">2</span>*np.random.rand(m,<span style="color: #B452CD">1</span>)
y = <span style="color: #B452CD">4</span>+<span style="color: #B452CD">3</span>*x+np.random.randn(m,<span style="color: #B452CD">1</span>)
xb = np.c_[np.ones((<span style="color: #B452CD">100</span>,<span style="color: #B452CD">1</span>)), x]
xb = np.c_[np.ones((m,<span style="color: #B452CD">1</span>)), x]
beta_linreg = np.linalg.inv(xb.T.dot(xb)).dot(xb.T).dot(y)
<span style="color: #8B008B; font-weight: bold">print</span>(beta_linreg)
beta = np.random.randn(<span style="color: #B452CD">2</span>,<span style="color: #B452CD">1</span>)
eta = <span style="color: #B452CD">0.1</span>
Niterations = <span style="color: #B452CD">1000</span>
m = <span style="color: #B452CD">100</span>
<span style="color: #8B008B; font-weight: bold">for</span> <span style="color: #658b00">iter</span> <span style="color: #8B008B">in</span> <span style="color: #658b00">range</span>(Niterations):
gradients = <span style="color: #B452CD">2.0</span>/m*xb.T.dot(xb.dot(beta)-y)
@@ -868,43 +848,11 @@ $$
$$
<p>
We can now extend our program to minimize \( C_{\text{ridge}}(\beta) \) using gradient descent and compare with the analytical solution given by
We can easily extend our program to minimize \( C_{\text{ridge}}(\beta) \) using gradient descent and compare with the analytical solution given by
$$
\beta_{\text{ridge}} = \left(X^T X + \lambda I_{2 \times 2} \right)^{-1} X^T \mathbf{y},
\beta_{\text{ridge}} = \left(X^T X + \lambda I_{2 \times 2} \right)^{-1} X^T \mathbf{y}.
$$
for \( \lambda = {0,1,10,50,100} \) (\( \lambda = 0 \) corresponds to ordinary least squares).
We can then compute \( ||\beta_{\text{ridge}}|| \) for each \( \lambda \).
<p>
<!-- code=python (!bc pycod) typeset with pygments style "perldoc" -->
<div class="highlight" style="background: #eeeedd"><pre style="line-height: 125%"><span></span><span style="color: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">numpy</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">np</span>
<span style="color: #CD5555">&quot;&quot;&quot;</span>
<span style="color: #CD5555">The following setup is just a suggestion, feel free to write it the way you like.</span>
<span style="color: #CD5555">&quot;&quot;&quot;</span>
<span style="color: #228B22">#Setup problem described in the exercise</span>
N = <span style="color: #B452CD">100</span> <span style="color: #228B22">#Nr of datapoints</span>
M = <span style="color: #B452CD">2</span> <span style="color: #228B22">#Nr of features</span>
x = np.random.rand(N)
y = <span style="color: #B452CD">5</span>*x**<span style="color: #B452CD">2</span> + <span style="color: #B452CD">0.1</span>*np.random.randn(N)
<span style="color: #228B22">#Compute analytic beta for Ridge regression </span>
X = np.c_[np.ones(N),x]
XT_X = np.dot(X.T,X)
l = <span style="color: #B452CD">0.1</span> <span style="color: #228B22">#Ridge parameter lambda</span>
Id = np.eye(XT_X.shape[<span style="color: #B452CD">0</span>])
Z = np.linalg.inv(XT_X+l*Id)
beta_ridge = np.dot(Z,np.dot(X.T,y))
<span style="color: #8B008B; font-weight: bold">print</span>(beta_ridge)
<span style="color: #8B008B; font-weight: bold">print</span>(np.linalg.norm(beta_ridge)) <span style="color: #228B22">#||beta||</span>
</pre></div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
+8 -60
View File
@@ -215,7 +215,7 @@ MathJax.Hub.Config({
<center>[2] <b>Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University</b></center>
<br>
<p>
<center><h4>Sep 19, 2019</h4></center> <!-- date -->
<center><h4>Sep 20, 2019</h4></center> <!-- date -->
<br>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
@@ -758,28 +758,7 @@ when \( ||\nabla_\beta C(\beta_k) || \leq \epsilon = 10^{-8} \).
<p>
And finally we can compare our solution for \( \beta \) with the analytic result given by
\( \beta= (X^TX)^{-1} X^T \mathbf{y} \).
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span><span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">numpy</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">np</span>
<span style="color: #BA2121; font-style: italic">&quot;&quot;&quot;</span>
<span style="color: #BA2121; font-style: italic">The following setup is just a suggestion, feel free to write it the way you like.</span>
<span style="color: #BA2121; font-style: italic">&quot;&quot;&quot;</span>
<span style="color: #408080; font-style: italic">#Setup problem described in the exercise</span>
N <span style="color: #666666">=</span> <span style="color: #666666">100</span> <span style="color: #408080; font-style: italic">#Nr of datapoints</span>
M <span style="color: #666666">=</span> <span style="color: #666666">2</span> <span style="color: #408080; font-style: italic">#Nr of features</span>
x <span style="color: #666666">=</span> np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>rand(N) <span style="color: #408080; font-style: italic">#Uniformly generated x-values in [0,1]</span>
y <span style="color: #666666">=</span> <span style="color: #666666">5*</span>x<span style="color: #666666">**2</span> <span style="color: #666666">+</span> <span style="color: #666666">0.1*</span>np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>randn(N)
X <span style="color: #666666">=</span> np<span style="color: #666666">.</span>c_[np<span style="color: #666666">.</span>ones(N),x] <span style="color: #408080; font-style: italic">#Construct design matrix</span>
<span style="color: #408080; font-style: italic">#Compute beta according to normal equations to compare with GD solution</span>
Xt_X_inv <span style="color: #666666">=</span> np<span style="color: #666666">.</span>linalg<span style="color: #666666">.</span>inv(np<span style="color: #666666">.</span>dot(X<span style="color: #666666">.</span>T,X))
Xt_y <span style="color: #666666">=</span> np<span style="color: #666666">.</span>dot(X<span style="color: #666666">.</span>transpose(),y)
beta_NE <span style="color: #666666">=</span> np<span style="color: #666666">.</span>dot(Xt_X_inv,Xt_y)
<span style="color: #008000; font-weight: bold">print</span>(beta_NE)
</pre></div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
@@ -799,17 +778,18 @@ Another simple example is here
<span style="color: #008000; font-weight: bold">from</span> <span style="color: #0000FF; font-weight: bold">matplotlib.ticker</span> <span style="color: #008000; font-weight: bold">import</span> LinearLocator, FormatStrFormatter
<span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">sys</span>
x <span style="color: #666666">=</span> <span style="color: #666666">2*</span>np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>rand(<span style="color: #666666">100</span>,<span style="color: #666666">1</span>)
y <span style="color: #666666">=</span> <span style="color: #666666">4+3*</span>x<span style="color: #666666">+</span>np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>randn(<span style="color: #666666">100</span>,<span style="color: #666666">1</span>)
<span style="color: #408080; font-style: italic"># the number of datapoints</span>
m <span style="color: #666666">=</span> <span style="color: #666666">100</span>
x <span style="color: #666666">=</span> <span style="color: #666666">2*</span>np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>rand(m,<span style="color: #666666">1</span>)
y <span style="color: #666666">=</span> <span style="color: #666666">4+3*</span>x<span style="color: #666666">+</span>np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>randn(m,<span style="color: #666666">1</span>)
xb <span style="color: #666666">=</span> np<span style="color: #666666">.</span>c_[np<span style="color: #666666">.</span>ones((<span style="color: #666666">100</span>,<span style="color: #666666">1</span>)), x]
xb <span style="color: #666666">=</span> np<span style="color: #666666">.</span>c_[np<span style="color: #666666">.</span>ones((m,<span style="color: #666666">1</span>)), x]
beta_linreg <span style="color: #666666">=</span> np<span style="color: #666666">.</span>linalg<span style="color: #666666">.</span>inv(xb<span style="color: #666666">.</span>T<span style="color: #666666">.</span>dot(xb))<span style="color: #666666">.</span>dot(xb<span style="color: #666666">.</span>T)<span style="color: #666666">.</span>dot(y)
<span style="color: #008000; font-weight: bold">print</span>(beta_linreg)
beta <span style="color: #666666">=</span> np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>randn(<span style="color: #666666">2</span>,<span style="color: #666666">1</span>)
eta <span style="color: #666666">=</span> <span style="color: #666666">0.1</span>
Niterations <span style="color: #666666">=</span> <span style="color: #666666">1000</span>
m <span style="color: #666666">=</span> <span style="color: #666666">100</span>
<span style="color: #008000; font-weight: bold">for</span> <span style="color: #008000">iter</span> <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #008000">range</span>(Niterations):
gradients <span style="color: #666666">=</span> <span style="color: #666666">2.0/</span>m<span style="color: #666666">*</span>xb<span style="color: #666666">.</span>T<span style="color: #666666">.</span>dot(xb<span style="color: #666666">.</span>dot(beta)<span style="color: #666666">-</span>y)
@@ -873,43 +853,11 @@ $$
$$
<p>
We can now extend our program to minimize \( C_{\text{ridge}}(\beta) \) using gradient descent and compare with the analytical solution given by
We can easily extend our program to minimize \( C_{\text{ridge}}(\beta) \) using gradient descent and compare with the analytical solution given by
$$
\beta_{\text{ridge}} = \left(X^T X + \lambda I_{2 \times 2} \right)^{-1} X^T \mathbf{y},
\beta_{\text{ridge}} = \left(X^T X + \lambda I_{2 \times 2} \right)^{-1} X^T \mathbf{y}.
$$
for \( \lambda = {0,1,10,50,100} \) (\( \lambda = 0 \) corresponds to ordinary least squares).
We can then compute \( ||\beta_{\text{ridge}}|| \) for each \( \lambda \).
<p>
<!-- code=python (!bc pycod) typeset with pygments style "default" -->
<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span><span style="color: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">numpy</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">np</span>
<span style="color: #BA2121; font-style: italic">&quot;&quot;&quot;</span>
<span style="color: #BA2121; font-style: italic">The following setup is just a suggestion, feel free to write it the way you like.</span>
<span style="color: #BA2121; font-style: italic">&quot;&quot;&quot;</span>
<span style="color: #408080; font-style: italic">#Setup problem described in the exercise</span>
N <span style="color: #666666">=</span> <span style="color: #666666">100</span> <span style="color: #408080; font-style: italic">#Nr of datapoints</span>
M <span style="color: #666666">=</span> <span style="color: #666666">2</span> <span style="color: #408080; font-style: italic">#Nr of features</span>
x <span style="color: #666666">=</span> np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>rand(N)
y <span style="color: #666666">=</span> <span style="color: #666666">5*</span>x<span style="color: #666666">**2</span> <span style="color: #666666">+</span> <span style="color: #666666">0.1*</span>np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>randn(N)
<span style="color: #408080; font-style: italic">#Compute analytic beta for Ridge regression </span>
X <span style="color: #666666">=</span> np<span style="color: #666666">.</span>c_[np<span style="color: #666666">.</span>ones(N),x]
XT_X <span style="color: #666666">=</span> np<span style="color: #666666">.</span>dot(X<span style="color: #666666">.</span>T,X)
l <span style="color: #666666">=</span> <span style="color: #666666">0.1</span> <span style="color: #408080; font-style: italic">#Ridge parameter lambda</span>
Id <span style="color: #666666">=</span> np<span style="color: #666666">.</span>eye(XT_X<span style="color: #666666">.</span>shape[<span style="color: #666666">0</span>])
Z <span style="color: #666666">=</span> np<span style="color: #666666">.</span>linalg<span style="color: #666666">.</span>inv(XT_X<span style="color: #666666">+</span>l<span style="color: #666666">*</span>Id)
beta_ridge <span style="color: #666666">=</span> np<span style="color: #666666">.</span>dot(Z,np<span style="color: #666666">.</span>dot(X<span style="color: #666666">.</span>T,y))
<span style="color: #008000; font-weight: bold">print</span>(beta_ridge)
<span style="color: #008000; font-weight: bold">print</span>(np<span style="color: #666666">.</span>linalg<span style="color: #666666">.</span>norm(beta_ridge)) <span style="color: #408080; font-style: italic">#||beta||</span>
</pre></div>
<p>
<!-- !split --><br><br><br><br><br><br><br><br><br><br>
+32 -107
View File
@@ -10,7 +10,7 @@
"<!-- Author: --> \n",
"**Morten Hjorth-Jensen**, Department of Physics, University of Oslo and Department of Physics and Astronomy and National Superconducting Cyclotron Laboratory, Michigan State University\n",
"\n",
"Date: **Sep 19, 2019**\n",
"Date: **Sep 20, 2019**\n",
"\n",
"Copyright 1999-2019, Morten Hjorth-Jensen. Released under CC Attribution-NonCommercial 4.0 license\n",
"\n",
@@ -717,41 +717,8 @@
"when $||\\nabla_\\beta C(\\beta_k) || \\leq \\epsilon = 10^{-8}$. \n",
"\n",
"And finally we can compare our solution for $\\beta$ with the analytic result given by \n",
"$\\beta= (X^TX)^{-1} X^T \\mathbf{y}$."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import numpy as np\n",
"$\\beta= (X^TX)^{-1} X^T \\mathbf{y}$.\n",
"\n",
"\"\"\"\n",
"The following setup is just a suggestion, feel free to write it the way you like.\n",
"\"\"\"\n",
"\n",
"#Setup problem described in the exercise\n",
"N = 100 #Nr of datapoints\n",
"M = 2 #Nr of features\n",
"x = np.random.rand(N) #Uniformly generated x-values in [0,1]\n",
"y = 5*x**2 + 0.1*np.random.randn(N)\n",
"X = np.c_[np.ones(N),x] #Construct design matrix\n",
"\n",
"#Compute beta according to normal equations to compare with GD solution\n",
"Xt_X_inv = np.linalg.inv(np.dot(X.T,X))\n",
"Xt_y = np.dot(X.transpose(),y)\n",
"beta_NE = np.dot(Xt_X_inv,Xt_y)\n",
"print(beta_NE)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Gradient Descent Example\n",
"\n",
"Another simple example is here"
@@ -759,7 +726,7 @@
},
{
"cell_type": "code",
"execution_count": 2,
"execution_count": 1,
"metadata": {
"collapsed": false
},
@@ -777,17 +744,18 @@
"from matplotlib.ticker import LinearLocator, FormatStrFormatter\n",
"import sys\n",
"\n",
"x = 2*np.random.rand(100,1)\n",
"y = 4+3*x+np.random.randn(100,1)\n",
"# the number of datapoints\n",
"m = 100\n",
"x = 2*np.random.rand(m,1)\n",
"y = 4+3*x+np.random.randn(m,1)\n",
"\n",
"xb = np.c_[np.ones((100,1)), x]\n",
"xb = np.c_[np.ones((m,1)), x]\n",
"beta_linreg = np.linalg.inv(xb.T.dot(xb)).dot(xb.T).dot(y)\n",
"print(beta_linreg)\n",
"beta = np.random.randn(2,1)\n",
"\n",
"eta = 0.1\n",
"Niterations = 1000\n",
"m = 100\n",
"\n",
"for iter in range(Niterations):\n",
" gradients = 2.0/m*xb.T.dot(xb.dot(beta)-y)\n",
@@ -817,7 +785,7 @@
},
{
"cell_type": "code",
"execution_count": 3,
"execution_count": 2,
"metadata": {
"collapsed": false
},
@@ -881,7 +849,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"We can now extend our program to minimize $C_{\\text{ridge}}(\\beta)$ using gradient descent and compare with the analytical solution given by"
"We can easily extend our program to minimize $C_{\\text{ridge}}(\\beta)$ using gradient descent and compare with the analytical solution given by"
]
},
{
@@ -889,53 +857,10 @@
"metadata": {},
"source": [
"$$\n",
"\\beta_{\\text{ridge}} = \\left(X^T X + \\lambda I_{2 \\times 2} \\right)^{-1} X^T \\mathbf{y},\n",
"\\beta_{\\text{ridge}} = \\left(X^T X + \\lambda I_{2 \\times 2} \\right)^{-1} X^T \\mathbf{y}.\n",
"$$"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"for $\\lambda = {0,1,10,50,100}$ ($\\lambda = 0$ corresponds to ordinary least squares). \n",
"We can then compute $||\\beta_{\\text{ridge}}||$ for each $\\lambda$."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"import numpy as np\n",
"\n",
"\"\"\"\n",
"The following setup is just a suggestion, feel free to write it the way you like.\n",
"\"\"\"\n",
"\n",
"#Setup problem described in the exercise\n",
"N = 100 #Nr of datapoints\n",
"M = 2 #Nr of features\n",
"x = np.random.rand(N)\n",
"y = 5*x**2 + 0.1*np.random.randn(N)\n",
"\n",
"\n",
"#Compute analytic beta for Ridge regression \n",
"X = np.c_[np.ones(N),x]\n",
"XT_X = np.dot(X.T,X)\n",
"\n",
"l = 0.1 #Ridge parameter lambda\n",
"Id = np.eye(XT_X.shape[0])\n",
"\n",
"Z = np.linalg.inv(XT_X+l*Id)\n",
"beta_ridge = np.dot(Z,np.dot(X.T,y))\n",
"\n",
"print(beta_ridge)\n",
"print(np.linalg.norm(beta_ridge)) #||beta||"
]
},
{
"cell_type": "markdown",
"metadata": {},
@@ -1051,7 +976,7 @@
},
{
"cell_type": "code",
"execution_count": 5,
"execution_count": 3,
"metadata": {
"collapsed": false
},
@@ -1115,7 +1040,7 @@
},
{
"cell_type": "code",
"execution_count": 6,
"execution_count": 4,
"metadata": {
"collapsed": false
},
@@ -1156,7 +1081,7 @@
},
{
"cell_type": "code",
"execution_count": 7,
"execution_count": 5,
"metadata": {
"collapsed": false
},
@@ -1633,7 +1558,7 @@
},
{
"cell_type": "code",
"execution_count": 8,
"execution_count": 6,
"metadata": {
"collapsed": false
},
@@ -1691,7 +1616,7 @@
},
{
"cell_type": "code",
"execution_count": 9,
"execution_count": 7,
"metadata": {
"collapsed": false
},
@@ -1729,7 +1654,7 @@
},
{
"cell_type": "code",
"execution_count": 10,
"execution_count": 8,
"metadata": {
"collapsed": false
},
@@ -1782,7 +1707,7 @@
},
{
"cell_type": "code",
"execution_count": 11,
"execution_count": 9,
"metadata": {
"collapsed": false
},
@@ -1824,7 +1749,7 @@
},
{
"cell_type": "code",
"execution_count": 12,
"execution_count": 10,
"metadata": {
"collapsed": false
},
@@ -1858,7 +1783,7 @@
},
{
"cell_type": "code",
"execution_count": 13,
"execution_count": 11,
"metadata": {
"collapsed": false
},
@@ -1892,7 +1817,7 @@
"metadata": {},
"source": [
"1\n",
"3\n",
"1\n",
" \n",
"<\n",
"<\n",
@@ -1920,7 +1845,7 @@
},
{
"cell_type": "code",
"execution_count": 14,
"execution_count": 12,
"metadata": {
"collapsed": false
},
@@ -1946,7 +1871,7 @@
},
{
"cell_type": "code",
"execution_count": 15,
"execution_count": 13,
"metadata": {
"collapsed": false
},
@@ -1995,7 +1920,7 @@
},
{
"cell_type": "code",
"execution_count": 16,
"execution_count": 14,
"metadata": {
"collapsed": false
},
@@ -2025,7 +1950,7 @@
},
{
"cell_type": "code",
"execution_count": 17,
"execution_count": 15,
"metadata": {
"collapsed": false
},
@@ -2055,7 +1980,7 @@
},
{
"cell_type": "code",
"execution_count": 18,
"execution_count": 16,
"metadata": {
"collapsed": false
},
@@ -2087,7 +2012,7 @@
},
{
"cell_type": "code",
"execution_count": 19,
"execution_count": 17,
"metadata": {
"collapsed": false
},
@@ -2410,7 +2335,7 @@
},
{
"cell_type": "code",
"execution_count": 20,
"execution_count": 18,
"metadata": {
"collapsed": false
},
@@ -2447,7 +2372,7 @@
},
{
"cell_type": "code",
"execution_count": 21,
"execution_count": 19,
"metadata": {
"collapsed": false
},
@@ -2467,7 +2392,7 @@
},
{
"cell_type": "code",
"execution_count": 22,
"execution_count": 20,
"metadata": {
"collapsed": false
},
@@ -2486,7 +2411,7 @@
},
{
"cell_type": "code",
"execution_count": 23,
"execution_count": 21,
"metadata": {
"collapsed": false
},
@@ -2510,7 +2435,7 @@
},
{
"cell_type": "code",
"execution_count": 24,
"execution_count": 22,
"metadata": {
"collapsed": false
},
Binary file not shown.
Binary file not shown.
+7 -55
View File
@@ -478,26 +478,6 @@ when $||\nabla_\beta C(\beta_k) || \leq \epsilon = 10^{-8}$.
And finally we can compare our solution for $\beta$ with the analytic result given by
$\beta= (X^TX)^{-1} X^T \mathbf{y}$.
!bc pycod
import numpy as np
"""
The following setup is just a suggestion, feel free to write it the way you like.
"""
#Setup problem described in the exercise
N = 100 #Nr of datapoints
M = 2 #Nr of features
x = np.random.rand(N) #Uniformly generated x-values in [0,1]
y = 5*x**2 + 0.1*np.random.randn(N)
X = np.c_[np.ones(N),x] #Construct design matrix
#Compute beta according to normal equations to compare with GD solution
Xt_X_inv = np.linalg.inv(np.dot(X.T,X))
Xt_y = np.dot(X.transpose(),y)
beta_NE = np.dot(Xt_X_inv,Xt_y)
print(beta_NE)
!ec
!split
===== Gradient Descent Example =====
@@ -514,17 +494,18 @@ from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import sys
x = 2*np.random.rand(100,1)
y = 4+3*x+np.random.randn(100,1)
# the number of datapoints
m = 100
x = 2*np.random.rand(m,1)
y = 4+3*x+np.random.randn(m,1)
xb = np.c_[np.ones((100,1)), x]
xb = np.c_[np.ones((m,1)), x]
beta_linreg = np.linalg.inv(xb.T.dot(xb)).dot(xb.T).dot(y)
print(beta_linreg)
beta = np.random.randn(2,1)
eta = 0.1
Niterations = 1000
m = 100
for iter in range(Niterations):
gradients = 2.0/m*xb.T.dot(xb.dot(beta)-y)
@@ -589,42 +570,13 @@ In order to minimize $C_{\text{ridge}}(\beta)$ using GD we only have adjust the
\]
!et
We can now extend our program to minimize $C_{\text{ridge}}(\beta)$ using gradient descent and compare with the analytical solution given by
We can easily extend our program to minimize $C_{\text{ridge}}(\beta)$ using gradient descent and compare with the analytical solution given by
!bt
\[
\beta_{\text{ridge}} = \left(X^T X + \lambda I_{2 \times 2} \right)^{-1} X^T \mathbf{y},
\beta_{\text{ridge}} = \left(X^T X + \lambda I_{2 \times 2} \right)^{-1} X^T \mathbf{y}.
\]
!et
for $\lambda = {0,1,10,50,100}$ ($\lambda = 0$ corresponds to ordinary least squares).
We can then compute $||\beta_{\text{ridge}}||$ for each $\lambda$.
!bc pycod
import numpy as np
"""
The following setup is just a suggestion, feel free to write it the way you like.
"""
#Setup problem described in the exercise
N = 100 #Nr of datapoints
M = 2 #Nr of features
x = np.random.rand(N)
y = 5*x**2 + 0.1*np.random.randn(N)
#Compute analytic beta for Ridge regression
X = np.c_[np.ones(N),x]
XT_X = np.dot(X.T,X)
l = 0.1 #Ridge parameter lambda
Id = np.eye(XT_X.shape[0])
Z = np.linalg.inv(XT_X+l*Id)
beta_ridge = np.dot(Z,np.dot(X.T,y))
print(beta_ridge)
print(np.linalg.norm(beta_ridge)) #||beta||
!ec
!split
===== Stochastic Gradient Descent =====