From ac368dd0c1a913ae43eaf7910113c3d1f3bd0d4e Mon Sep 17 00:00:00 2001
From: mhjensen <morten.hjorth-jensen@fys.uio.no>
Date: Thu, 2 Jan 2020 10:45:16 +0100
Subject: [PATCH] more pca

---
 doc/pub/DimRed/html/._DimRed-bs018.html      |  41 +++++++++++-------
 doc/pub/DimRed/html/DimRed-reveal.html       |  40 +++++++++++-------
 doc/pub/DimRed/html/DimRed-solarized.html    |  41 +++++++++++-------
 doc/pub/DimRed/html/DimRed.html              |  41 +++++++++++-------
 doc/pub/DimRed/ipynb/DimRed.ipynb            |  40 +++++++++++-------
 doc/pub/DimRed/ipynb/ipynb-DimRed-src.tar.gz | Bin 190 -> 191 bytes
 doc/pub/DimRed/pdf/DimRed-minted.pdf         | Bin 263343 -> 263641 bytes
 doc/src/DimRed/DimRed.do.txt                 |  42 +++++++++++--------
 doc/src/DimRed/Programs/PCAsimple.py         |  41 +++++++++---------
 9 files changed, 177 insertions(+), 109 deletions(-)

diff --git a/doc/pub/DimRed/html/._DimRed-bs018.html b/doc/pub/DimRed/html/._DimRed-bs018.html
index 26757403b..a7cce97b1 100644
--- a/doc/pub/DimRed/html/._DimRed-bs018.html
+++ b/doc/pub/DimRed/html/._DimRed-bs018.html
@@ -229,7 +229,7 @@ Note that the function <b>multivariate</b> returns also the covariance discussed
 <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: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">pandas</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">pd</span>
 <span style="color: #008000; font-weight: bold">from</span> <span style="color: #0000FF; font-weight: bold">IPython.display</span> <span style="color: #008000; font-weight: bold">import</span> display
-n <span style="color: #666666">=</span> <span style="color: #666666">100</span>
+n <span style="color: #666666">=</span> <span style="color: #666666">10000</span>
 mean <span style="color: #666666">=</span> (<span style="color: #666666">-1</span>, <span style="color: #666666">2</span>)
 cov <span style="color: #666666">=</span> [[<span style="color: #666666">4</span>, <span style="color: #666666">2</span>], [<span style="color: #666666">2</span>, <span style="color: #666666">2</span>]]
 X <span style="color: #666666">=</span> np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>multivariate_normal(mean, cov, n)
@@ -313,14 +313,15 @@ plt<span style="color: #666666">.</span>show()
 </pre></div>
 <p>
 Depending on the number of points \( n \), we will get results that are close to the covariance values defined above.
+The plot shows how the data are clustered around a line with slope close to one. Is this expected?
 
 <h3 id="___sec20" class="anchor">Diagonalize the sample covariance matrix to obtain the principal components </h3>
 
 <p>
 Now we are ready to solve for the principal components! To do so we
-diagonalize the sample covariance matrix \( \Sigma_n \). We can use the
+diagonalize the sample covariance matrix \( \Sigma \). We can use the
 function <b>np.linalg.eig</b> to do so. It will return the eigenvalues and
-eigenvectors of \( \Sigma_n \). Once we have these we can perform the 
+eigenvectors of \( \Sigma \). Once we have these we can perform the 
 following tasks:
 
 <ul>
@@ -343,26 +344,38 @@ Collecting all these steps we can write our own PCA function and
 compare this with the functionality included in <b>Scikit-Learn</b>.
 
 <p>
-The code here outlines some of the elements we could include in the analysis. Feel free to extend upon this.
+The code here outlines some of the elements we could include in the
+analysis. Feel free to extend upon this in order to address the above
+questions.
+
 <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: #408080; font-style: italic">#Now we do an SVD</span>
-U, s, V <span style="color: #666666">=</span> np<span style="color: #666666">.</span>linalg<span style="color: #666666">.</span>svd(X_centered)
-c1 <span style="color: #666666">=</span> V<span style="color: #666666">.</span>T[:, <span style="color: #666666">0</span>]
-c2 <span style="color: #666666">=</span> V<span style="color: #666666">.</span>T[:, <span style="color: #666666">1</span>]
-W2 <span style="color: #666666">=</span> V<span style="color: #666666">.</span>T[:, :<span style="color: #666666">2</span>]
-X2D <span style="color: #666666">=</span> X_centered<span style="color: #666666">.</span>dot(W2)
-<span style="color: #008000; font-weight: bold">print</span>(X2D)
+<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span><span style="color: #408080; font-style: italic"># diagonalize and obtain eigenvalues, not necessarily sorted</span>
+EigValues, EigVectors <span style="color: #666666">=</span> np<span style="color: #666666">.</span>linalg<span style="color: #666666">.</span>eig(Cov)
+<span style="color: #408080; font-style: italic"># sort eigenvectors and eigenvalues</span>
+<span style="color: #408080; font-style: italic">#permute = EigValues.argsort()</span>
+<span style="color: #408080; font-style: italic">#EigValues = EigValues[permute]</span>
+<span style="color: #408080; font-style: italic">#EigVectors = EigVectors[:,permute]</span>
+<span style="color: #008000; font-weight: bold">print</span>(<span style="color: #BA2121">&quot;Eigenvalues of Covariance matrix&quot;</span>)
+<span style="color: #008000; font-weight: bold">for</span> i <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #008000">range</span>(<span style="color: #666666">2</span>):
+    <span style="color: #008000; font-weight: bold">print</span>(EigValues[i])
+FirstEigvector <span style="color: #666666">=</span> EigVectors[:,<span style="color: #666666">0</span>]
+SecondEigvector <span style="color: #666666">=</span> EigVectors[:,<span style="color: #666666">1</span>]
+<span style="color: #008000; font-weight: bold">print</span>(<span style="color: #BA2121">&quot;First eigenvector&quot;</span>)
+<span style="color: #008000; font-weight: bold">print</span>(FirstEigvector)
+<span style="color: #008000; font-weight: bold">print</span>(<span style="color: #BA2121">&quot;Second eigenvector&quot;</span>)
+<span style="color: #008000; font-weight: bold">print</span>(SecondEigvector)
 <span style="color: #408080; font-style: italic">#thereafter we do a PCA with Scikit-learn</span>
 <span style="color: #008000; font-weight: bold">from</span> <span style="color: #0000FF; font-weight: bold">sklearn.decomposition</span> <span style="color: #008000; font-weight: bold">import</span> PCA
 pca <span style="color: #666666">=</span> PCA(n_components <span style="color: #666666">=</span> <span style="color: #666666">2</span>)
 X2Dsl <span style="color: #666666">=</span> pca<span style="color: #666666">.</span>fit_transform(X)
-<span style="color: #008000; font-weight: bold">print</span>(<span style="color: #BA2121">&quot;Check that we get the same&quot;</span>)
-<span style="color: #008000; font-weight: bold">print</span>(X2D<span style="color: #666666">-</span>X2Dsl)
-
+<span style="color: #008000; font-weight: bold">print</span>(<span style="color: #BA2121">&quot;Eigenvector of largest eigenvalue&quot;</span>)
 <span style="color: #008000; font-weight: bold">print</span>(pca<span style="color: #666666">.</span>components_<span style="color: #666666">.</span>T[:, <span style="color: #666666">0</span>])
 </pre></div>
+<p>
+This code does not contain all the above elements, but it shows how we can use <b>Scikit-Learn</b> to extract the eigenvector which corresponds to the largest eigenvalue. Try to, based on the above, to address the questions above.
+
 <p>
 <p>
 <!-- navigation buttons at the bottom of the page -->
diff --git a/doc/pub/DimRed/html/DimRed-reveal.html b/doc/pub/DimRed/html/DimRed-reveal.html
index 5a9a90a8f..367451ae8 100644
--- a/doc/pub/DimRed/html/DimRed-reveal.html
+++ b/doc/pub/DimRed/html/DimRed-reveal.html
@@ -1035,7 +1035,7 @@ Note that the function <b>multivariate</b> returns also the covariance discussed
 <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: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">pandas</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">pd</span>
 <span style="color: #8B008B; font-weight: bold">from</span> <span style="color: #008b45; text-decoration: underline">IPython.display</span> <span style="color: #8B008B; font-weight: bold">import</span> display
-n = <span style="color: #B452CD">100</span>
+n = <span style="color: #B452CD">10000</span>
 mean = (-<span style="color: #B452CD">1</span>, <span style="color: #B452CD">2</span>)
 cov = [[<span style="color: #B452CD">4</span>, <span style="color: #B452CD">2</span>], [<span style="color: #B452CD">2</span>, <span style="color: #B452CD">2</span>]]
 X = np.random.multivariate_normal(mean, cov, n)
@@ -1125,14 +1125,15 @@ plt.show()
 </pre></div>
 <p>
 Depending on the number of points \( n \), we will get results that are close to the covariance values defined above.
+The plot shows how the data are clustered around a line with slope close to one. Is this expected?
 
 <h3 id="___sec20">Diagonalize the sample covariance matrix to obtain the principal components </h3>
 
 <p>
 Now we are ready to solve for the principal components! To do so we
-diagonalize the sample covariance matrix \( \Sigma_n \). We can use the
+diagonalize the sample covariance matrix \( \Sigma \). We can use the
 function <b>np.linalg.eig</b> to do so. It will return the eigenvalues and
-eigenvectors of \( \Sigma_n \). Once we have these we can perform the 
+eigenvectors of \( \Sigma \). Once we have these we can perform the 
 following tasks:
 
 <ul>
@@ -1156,26 +1157,37 @@ Collecting all these steps we can write our own PCA function and
 compare this with the functionality included in <b>Scikit-Learn</b>.
 
 <p>
-The code here outlines some of the elements we could include in the analysis. Feel free to extend upon this.
+The code here outlines some of the elements we could include in the
+analysis. Feel free to extend upon this in order to address the above
+questions.
+
 <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: #228B22">#Now we do an SVD</span>
-U, s, V = np.linalg.svd(X_centered)
-c1 = V.T[:, <span style="color: #B452CD">0</span>]
-c2 = V.T[:, <span style="color: #B452CD">1</span>]
-W2 = V.T[:, :<span style="color: #B452CD">2</span>]
-X2D = X_centered.dot(W2)
-<span style="color: #8B008B; font-weight: bold">print</span>(X2D)
+<div class="highlight" style="background: #eeeedd"><pre style="font-size: 80%; line-height: 125%"><span></span><span style="color: #228B22"># diagonalize and obtain eigenvalues, not necessarily sorted</span>
+EigValues, EigVectors = np.linalg.eig(Cov)
+<span style="color: #228B22"># sort eigenvectors and eigenvalues</span>
+<span style="color: #228B22">#permute = EigValues.argsort()</span>
+<span style="color: #228B22">#EigValues = EigValues[permute]</span>
+<span style="color: #228B22">#EigVectors = EigVectors[:,permute]</span>
+<span style="color: #8B008B; font-weight: bold">print</span>(<span style="color: #CD5555">&quot;Eigenvalues of Covariance matrix&quot;</span>)
+<span style="color: #8B008B; font-weight: bold">for</span> i <span style="color: #8B008B">in</span> <span style="color: #658b00">range</span>(<span style="color: #B452CD">2</span>):
+    <span style="color: #8B008B; font-weight: bold">print</span>(EigValues[i])
+FirstEigvector = EigVectors[:,<span style="color: #B452CD">0</span>]
+SecondEigvector = EigVectors[:,<span style="color: #B452CD">1</span>]
+<span style="color: #8B008B; font-weight: bold">print</span>(<span style="color: #CD5555">&quot;First eigenvector&quot;</span>)
+<span style="color: #8B008B; font-weight: bold">print</span>(FirstEigvector)
+<span style="color: #8B008B; font-weight: bold">print</span>(<span style="color: #CD5555">&quot;Second eigenvector&quot;</span>)
+<span style="color: #8B008B; font-weight: bold">print</span>(SecondEigvector)
 <span style="color: #228B22">#thereafter we do a PCA with Scikit-learn</span>
 <span style="color: #8B008B; font-weight: bold">from</span> <span style="color: #008b45; text-decoration: underline">sklearn.decomposition</span> <span style="color: #8B008B; font-weight: bold">import</span> PCA
 pca = PCA(n_components = <span style="color: #B452CD">2</span>)
 X2Dsl = pca.fit_transform(X)
-<span style="color: #8B008B; font-weight: bold">print</span>(<span style="color: #CD5555">&quot;Check that we get the same&quot;</span>)
-<span style="color: #8B008B; font-weight: bold">print</span>(X2D-X2Dsl)
-
+<span style="color: #8B008B; font-weight: bold">print</span>(<span style="color: #CD5555">&quot;Eigenvector of largest eigenvalue&quot;</span>)
 <span style="color: #8B008B; font-weight: bold">print</span>(pca.components_.T[:, <span style="color: #B452CD">0</span>])
 </pre></div>
+<p>
+This code does not contain all the above elements, but it shows how we can use <b>Scikit-Learn</b> to extract the eigenvector which corresponds to the largest eigenvalue. Try to, based on the above, to address the questions above.
 </section>
 
 
diff --git a/doc/pub/DimRed/html/DimRed-solarized.html b/doc/pub/DimRed/html/DimRed-solarized.html
index 0c1ae2da3..29d7191d6 100644
--- a/doc/pub/DimRed/html/DimRed-solarized.html
+++ b/doc/pub/DimRed/html/DimRed-solarized.html
@@ -1015,7 +1015,7 @@ Note that the function <b>multivariate</b> returns also the covariance discussed
 <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: #8B008B; font-weight: bold">import</span> <span style="color: #008b45; text-decoration: underline">pandas</span> <span style="color: #8B008B; font-weight: bold">as</span> <span style="color: #008b45; text-decoration: underline">pd</span>
 <span style="color: #8B008B; font-weight: bold">from</span> <span style="color: #008b45; text-decoration: underline">IPython.display</span> <span style="color: #8B008B; font-weight: bold">import</span> display
-n = <span style="color: #B452CD">100</span>
+n = <span style="color: #B452CD">10000</span>
 mean = (-<span style="color: #B452CD">1</span>, <span style="color: #B452CD">2</span>)
 cov = [[<span style="color: #B452CD">4</span>, <span style="color: #B452CD">2</span>], [<span style="color: #B452CD">2</span>, <span style="color: #B452CD">2</span>]]
 X = np.random.multivariate_normal(mean, cov, n)
@@ -1099,14 +1099,15 @@ plt.show()
 </pre></div>
 <p>
 Depending on the number of points \( n \), we will get results that are close to the covariance values defined above.
+The plot shows how the data are clustered around a line with slope close to one. Is this expected?
 
 <h3 id="___sec20">Diagonalize the sample covariance matrix to obtain the principal components </h3>
 
 <p>
 Now we are ready to solve for the principal components! To do so we
-diagonalize the sample covariance matrix \( \Sigma_n \). We can use the
+diagonalize the sample covariance matrix \( \Sigma \). We can use the
 function <b>np.linalg.eig</b> to do so. It will return the eigenvalues and
-eigenvectors of \( \Sigma_n \). Once we have these we can perform the 
+eigenvectors of \( \Sigma \). Once we have these we can perform the 
 following tasks:
 
 <ul>
@@ -1129,26 +1130,38 @@ Collecting all these steps we can write our own PCA function and
 compare this with the functionality included in <b>Scikit-Learn</b>.
 
 <p>
-The code here outlines some of the elements we could include in the analysis. Feel free to extend upon this.
+The code here outlines some of the elements we could include in the
+analysis. Feel free to extend upon this in order to address the above
+questions.
+
 <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: #228B22">#Now we do an SVD</span>
-U, s, V = np.linalg.svd(X_centered)
-c1 = V.T[:, <span style="color: #B452CD">0</span>]
-c2 = V.T[:, <span style="color: #B452CD">1</span>]
-W2 = V.T[:, :<span style="color: #B452CD">2</span>]
-X2D = X_centered.dot(W2)
-<span style="color: #8B008B; font-weight: bold">print</span>(X2D)
+<div class="highlight" style="background: #eeeedd"><pre style="line-height: 125%"><span></span><span style="color: #228B22"># diagonalize and obtain eigenvalues, not necessarily sorted</span>
+EigValues, EigVectors = np.linalg.eig(Cov)
+<span style="color: #228B22"># sort eigenvectors and eigenvalues</span>
+<span style="color: #228B22">#permute = EigValues.argsort()</span>
+<span style="color: #228B22">#EigValues = EigValues[permute]</span>
+<span style="color: #228B22">#EigVectors = EigVectors[:,permute]</span>
+<span style="color: #8B008B; font-weight: bold">print</span>(<span style="color: #CD5555">&quot;Eigenvalues of Covariance matrix&quot;</span>)
+<span style="color: #8B008B; font-weight: bold">for</span> i <span style="color: #8B008B">in</span> <span style="color: #658b00">range</span>(<span style="color: #B452CD">2</span>):
+    <span style="color: #8B008B; font-weight: bold">print</span>(EigValues[i])
+FirstEigvector = EigVectors[:,<span style="color: #B452CD">0</span>]
+SecondEigvector = EigVectors[:,<span style="color: #B452CD">1</span>]
+<span style="color: #8B008B; font-weight: bold">print</span>(<span style="color: #CD5555">&quot;First eigenvector&quot;</span>)
+<span style="color: #8B008B; font-weight: bold">print</span>(FirstEigvector)
+<span style="color: #8B008B; font-weight: bold">print</span>(<span style="color: #CD5555">&quot;Second eigenvector&quot;</span>)
+<span style="color: #8B008B; font-weight: bold">print</span>(SecondEigvector)
 <span style="color: #228B22">#thereafter we do a PCA with Scikit-learn</span>
 <span style="color: #8B008B; font-weight: bold">from</span> <span style="color: #008b45; text-decoration: underline">sklearn.decomposition</span> <span style="color: #8B008B; font-weight: bold">import</span> PCA
 pca = PCA(n_components = <span style="color: #B452CD">2</span>)
 X2Dsl = pca.fit_transform(X)
-<span style="color: #8B008B; font-weight: bold">print</span>(<span style="color: #CD5555">&quot;Check that we get the same&quot;</span>)
-<span style="color: #8B008B; font-weight: bold">print</span>(X2D-X2Dsl)
-
+<span style="color: #8B008B; font-weight: bold">print</span>(<span style="color: #CD5555">&quot;Eigenvector of largest eigenvalue&quot;</span>)
 <span style="color: #8B008B; font-weight: bold">print</span>(pca.components_.T[:, <span style="color: #B452CD">0</span>])
 </pre></div>
+<p>
+This code does not contain all the above elements, but it shows how we can use <b>Scikit-Learn</b> to extract the eigenvector which corresponds to the largest eigenvalue. Try to, based on the above, to address the questions above.
+
 <p>
 <!-- !split --><br><br><br><br><br><br><br><br><br><br>
 
diff --git a/doc/pub/DimRed/html/DimRed.html b/doc/pub/DimRed/html/DimRed.html
index 479e67b4f..8960e958f 100644
--- a/doc/pub/DimRed/html/DimRed.html
+++ b/doc/pub/DimRed/html/DimRed.html
@@ -1020,7 +1020,7 @@ Note that the function <b>multivariate</b> returns also the covariance discussed
 <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: #008000; font-weight: bold">import</span> <span style="color: #0000FF; font-weight: bold">pandas</span> <span style="color: #008000; font-weight: bold">as</span> <span style="color: #0000FF; font-weight: bold">pd</span>
 <span style="color: #008000; font-weight: bold">from</span> <span style="color: #0000FF; font-weight: bold">IPython.display</span> <span style="color: #008000; font-weight: bold">import</span> display
-n <span style="color: #666666">=</span> <span style="color: #666666">100</span>
+n <span style="color: #666666">=</span> <span style="color: #666666">10000</span>
 mean <span style="color: #666666">=</span> (<span style="color: #666666">-1</span>, <span style="color: #666666">2</span>)
 cov <span style="color: #666666">=</span> [[<span style="color: #666666">4</span>, <span style="color: #666666">2</span>], [<span style="color: #666666">2</span>, <span style="color: #666666">2</span>]]
 X <span style="color: #666666">=</span> np<span style="color: #666666">.</span>random<span style="color: #666666">.</span>multivariate_normal(mean, cov, n)
@@ -1104,14 +1104,15 @@ plt<span style="color: #666666">.</span>show()
 </pre></div>
 <p>
 Depending on the number of points \( n \), we will get results that are close to the covariance values defined above.
+The plot shows how the data are clustered around a line with slope close to one. Is this expected?
 
 <h3 id="___sec20">Diagonalize the sample covariance matrix to obtain the principal components </h3>
 
 <p>
 Now we are ready to solve for the principal components! To do so we
-diagonalize the sample covariance matrix \( \Sigma_n \). We can use the
+diagonalize the sample covariance matrix \( \Sigma \). We can use the
 function <b>np.linalg.eig</b> to do so. It will return the eigenvalues and
-eigenvectors of \( \Sigma_n \). Once we have these we can perform the 
+eigenvectors of \( \Sigma \). Once we have these we can perform the 
 following tasks:
 
 <ul>
@@ -1134,26 +1135,38 @@ Collecting all these steps we can write our own PCA function and
 compare this with the functionality included in <b>Scikit-Learn</b>.
 
 <p>
-The code here outlines some of the elements we could include in the analysis. Feel free to extend upon this.
+The code here outlines some of the elements we could include in the
+analysis. Feel free to extend upon this in order to address the above
+questions.
+
 <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: #408080; font-style: italic">#Now we do an SVD</span>
-U, s, V <span style="color: #666666">=</span> np<span style="color: #666666">.</span>linalg<span style="color: #666666">.</span>svd(X_centered)
-c1 <span style="color: #666666">=</span> V<span style="color: #666666">.</span>T[:, <span style="color: #666666">0</span>]
-c2 <span style="color: #666666">=</span> V<span style="color: #666666">.</span>T[:, <span style="color: #666666">1</span>]
-W2 <span style="color: #666666">=</span> V<span style="color: #666666">.</span>T[:, :<span style="color: #666666">2</span>]
-X2D <span style="color: #666666">=</span> X_centered<span style="color: #666666">.</span>dot(W2)
-<span style="color: #008000; font-weight: bold">print</span>(X2D)
+<div class="highlight" style="background: #f8f8f8"><pre style="line-height: 125%"><span></span><span style="color: #408080; font-style: italic"># diagonalize and obtain eigenvalues, not necessarily sorted</span>
+EigValues, EigVectors <span style="color: #666666">=</span> np<span style="color: #666666">.</span>linalg<span style="color: #666666">.</span>eig(Cov)
+<span style="color: #408080; font-style: italic"># sort eigenvectors and eigenvalues</span>
+<span style="color: #408080; font-style: italic">#permute = EigValues.argsort()</span>
+<span style="color: #408080; font-style: italic">#EigValues = EigValues[permute]</span>
+<span style="color: #408080; font-style: italic">#EigVectors = EigVectors[:,permute]</span>
+<span style="color: #008000; font-weight: bold">print</span>(<span style="color: #BA2121">&quot;Eigenvalues of Covariance matrix&quot;</span>)
+<span style="color: #008000; font-weight: bold">for</span> i <span style="color: #AA22FF; font-weight: bold">in</span> <span style="color: #008000">range</span>(<span style="color: #666666">2</span>):
+    <span style="color: #008000; font-weight: bold">print</span>(EigValues[i])
+FirstEigvector <span style="color: #666666">=</span> EigVectors[:,<span style="color: #666666">0</span>]
+SecondEigvector <span style="color: #666666">=</span> EigVectors[:,<span style="color: #666666">1</span>]
+<span style="color: #008000; font-weight: bold">print</span>(<span style="color: #BA2121">&quot;First eigenvector&quot;</span>)
+<span style="color: #008000; font-weight: bold">print</span>(FirstEigvector)
+<span style="color: #008000; font-weight: bold">print</span>(<span style="color: #BA2121">&quot;Second eigenvector&quot;</span>)
+<span style="color: #008000; font-weight: bold">print</span>(SecondEigvector)
 <span style="color: #408080; font-style: italic">#thereafter we do a PCA with Scikit-learn</span>
 <span style="color: #008000; font-weight: bold">from</span> <span style="color: #0000FF; font-weight: bold">sklearn.decomposition</span> <span style="color: #008000; font-weight: bold">import</span> PCA
 pca <span style="color: #666666">=</span> PCA(n_components <span style="color: #666666">=</span> <span style="color: #666666">2</span>)
 X2Dsl <span style="color: #666666">=</span> pca<span style="color: #666666">.</span>fit_transform(X)
-<span style="color: #008000; font-weight: bold">print</span>(<span style="color: #BA2121">&quot;Check that we get the same&quot;</span>)
-<span style="color: #008000; font-weight: bold">print</span>(X2D<span style="color: #666666">-</span>X2Dsl)
-
+<span style="color: #008000; font-weight: bold">print</span>(<span style="color: #BA2121">&quot;Eigenvector of largest eigenvalue&quot;</span>)
 <span style="color: #008000; font-weight: bold">print</span>(pca<span style="color: #666666">.</span>components_<span style="color: #666666">.</span>T[:, <span style="color: #666666">0</span>])
 </pre></div>
+<p>
+This code does not contain all the above elements, but it shows how we can use <b>Scikit-Learn</b> to extract the eigenvector which corresponds to the largest eigenvalue. Try to, based on the above, to address the questions above.
+
 <p>
 <!-- !split --><br><br><br><br><br><br><br><br><br><br>
 
diff --git a/doc/pub/DimRed/ipynb/DimRed.ipynb b/doc/pub/DimRed/ipynb/DimRed.ipynb
index bded76346..2301b9649 100644
--- a/doc/pub/DimRed/ipynb/DimRed.ipynb
+++ b/doc/pub/DimRed/ipynb/DimRed.ipynb
@@ -1108,7 +1108,7 @@
     "import numpy as np\n",
     "import pandas as pd\n",
     "from IPython.display import display\n",
-    "n = 100\n",
+    "n = 10000\n",
     "mean = (-1, 2)\n",
     "cov = [[4, 2], [2, 2]]\n",
     "X = np.random.multivariate_normal(mean, cov, n)"
@@ -1258,14 +1258,14 @@
    "metadata": {},
    "source": [
     "Depending on the number of points $n$, we will get results that are close to the covariance values defined above.\n",
-    "\n",
+    "The plot shows how the data are clustered around a line with slope close to one. Is this expected?\n",
     "\n",
     "### Diagonalize the sample covariance matrix to obtain the principal components\n",
     "\n",
     "Now we are ready to solve for the principal components! To do so we\n",
-    "diagonalize the sample covariance matrix $\\Sigma_n$. We can use the\n",
+    "diagonalize the sample covariance matrix $\\Sigma$. We can use the\n",
     "function **np.linalg.eig** to do so. It will return the eigenvalues and\n",
-    "eigenvectors of $\\Sigma_n$. Once we have these we can perform the \n",
+    "eigenvectors of $\\Sigma$. Once we have these we can perform the \n",
     "following tasks:\n",
     "\n",
     "* We compute the percentage of the total variance captured by the first principal component\n",
@@ -1295,7 +1295,9 @@
     "Collecting all these steps we can write our own PCA function and\n",
     "compare this with the functionality included in **Scikit-Learn**.  \n",
     "\n",
-    "The code here outlines some of the elements we could include in the analysis. Feel free to extend upon this."
+    "The code here outlines some of the elements we could include in the\n",
+    "analysis. Feel free to extend upon this in order to address the above\n",
+    "questions."
    ]
   },
   {
@@ -1306,20 +1308,26 @@
    },
    "outputs": [],
    "source": [
-    "#Now we do an SVD\n",
-    "U, s, V = np.linalg.svd(X_centered)\n",
-    "c1 = V.T[:, 0]\n",
-    "c2 = V.T[:, 1]\n",
-    "W2 = V.T[:, :2]\n",
-    "X2D = X_centered.dot(W2)\n",
-    "print(X2D)\n",
+    "# diagonalize and obtain eigenvalues, not necessarily sorted\n",
+    "EigValues, EigVectors = np.linalg.eig(Cov)\n",
+    "# sort eigenvectors and eigenvalues\n",
+    "#permute = EigValues.argsort()\n",
+    "#EigValues = EigValues[permute]\n",
+    "#EigVectors = EigVectors[:,permute]\n",
+    "print(\"Eigenvalues of Covariance matrix\")\n",
+    "for i in range(2):\n",
+    "    print(EigValues[i])\n",
+    "FirstEigvector = EigVectors[:,0]\n",
+    "SecondEigvector = EigVectors[:,1]\n",
+    "print(\"First eigenvector\")\n",
+    "print(FirstEigvector)\n",
+    "print(\"Second eigenvector\")\n",
+    "print(SecondEigvector)\n",
     "#thereafter we do a PCA with Scikit-learn\n",
     "from sklearn.decomposition import PCA\n",
     "pca = PCA(n_components = 2)\n",
     "X2Dsl = pca.fit_transform(X)\n",
-    "print(\"Check that we get the same\")\n",
-    "print(X2D-X2Dsl)\n",
-    "\n",
+    "print(\"Eigenvector of largest eigenvalue\")\n",
     "print(pca.components_.T[:, 0])"
    ]
   },
@@ -1327,6 +1335,8 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
+    "This code does not contain all the above elements, but it shows how we can use **Scikit-Learn** to extract the eigenvector which corresponds to the largest eigenvalue. Try to, based on the above, to address the questions above. \n",
+    "\n",
     "## Classical PCA Theorem\n",
     "\n",
     "We assume now that we have a design matrix $\\boldsymbol{X}$ which has been\n",
diff --git a/doc/pub/DimRed/ipynb/ipynb-DimRed-src.tar.gz b/doc/pub/DimRed/ipynb/ipynb-DimRed-src.tar.gz
index 603419a5eb62e8043df1425cde6a9db54ab251ac..88a78300f21e6f6ecd96e741e07dc05d57a8669c 100644
GIT binary patch
literal 191
zcmV;w06_mAiwFRsyA56d1MSaC3c@fD2H>uHia9|^nl`S5cHu&h;ssKYny5``l7hXx
zeSoeMH${Yen}0%vVWuj##lB4Z-AAiI2+=HqDK|8q6Q#?0hH(y<P(%pHX_f-QV?iQ7
z>#g+CJJ0O6iaZPLgnrvMwl&R%Im;>V%s+9I%E9h1cxg0f<3(-E1q#(xf+wmQbPe1H
tH+_K>w_aHROb@^~S6&&-u5oA0$+FfZ@vn~wf*=UK_5fC{sH^}8002ikRucdK

literal 190
zcmV;v073sBiwFRlwhdkY1MSbv3c@f92k@Qu6nTP?uG>5d+rfh%!x!ji=c=xqZ3ph{
z-3RDN@iIi{@A4-kgrqFC`M!wU-36;b#4Jf*%1xHeiSbg;2uXp^4CyGQ6p&DwqygZ1
zC%yF6F*BahoQ1MNy`3A!>iWZ;<rR45pEx$s!tOA*Mk`R-iAw7W<gzUVk7PHf3b+<_
sYJm~AUK#;R4Zt{;PHM%jVQ2Kou+~Q5ub=Te&-1?a0Az}Q695PR0IcF()Bpeg

diff --git a/doc/pub/DimRed/pdf/DimRed-minted.pdf b/doc/pub/DimRed/pdf/DimRed-minted.pdf
index 6c4542cb33528bc0b58deb310b029535bbf40484..c79cb258405889680f98e9f788a8e288193ed36f 100644
GIT binary patch
delta 5549
zcmV;e6;kT2i4fU~5U>qd1T!%+F_RHlDu3N?<hBui_g}Gn)1qFUD3TJ{A}HERu1HZN
zZDXJ<l04{kWm}C_+DO{-IS>8YJHsI*Qr6nW;ac=%m!kMFobUNj>E>>eZvJ?cy8l1i
zT)lb6i;ZC>&+^U9?S`{7<?t@u@C_G&6}jHzs$iK=o15e2`{Y=a?RF=GNIq;;o`00h
zxuVB<p#NX0zQb0Ql<k4K#``Lq(${6XU(D*-iYEP#s<isB&5ESj@}z&H*ZZzHAKkun
z8?Blphia(r+BdgVIiCA!`{T`@HYqcDc?zjmrrjaRv9CYER-A2jGBe3Pw^}9LdArN<
zJbBn=8Xg+h|A0M(-HInIJ+$3O<A3fd*RYgeV<?ZV)mOv08Ha0l(mWaJb`J|wv;eGc
z5;yKpf2a?x{T}yyTr3GJD?H2C+Sw%E<nSjlC!iZUQ>U<WdDpdNLr?n-cZ~&r{lp4t
zS50-Kg90JrK%EbD<CXxLT57$=zr?KoqQhgFEP1%En`#<Mjkard6D3Ig6MyZ>I1>fY
z4of+GXsgOAkJTn=_e7<z55zV%!to6jfAcOAp)zMF<2i14bAZ(|8X|ZwL;^GUM#C^o
zgHAZ;>Z_ZppRYWnATyS60hAh60IB=q)%QQ9n*+Q=X$TO<gB^Nw%cYxU^WD{7$-E+p
z=B4P+0a=k5Y9`Xm(|z0Z&VNYGJsp1lA)QDD`}&7e$h|X@av&vvZL?sx9haOHve@i0
z#k9$3-M8KGbRLnT+FO&(P=`;A>jb3%`}(lWWu$^DKEpU(i8N;lNLNNNSx6G=Hyujo
z0clK8`$>$Tj%8m_=UrE~cht1<hIi+|y>oUBs*1Wv-Usx#YIm#pq<^$viW}N_ta`5{
zhjJ`qjoMG=4}H{he(QFS*r#5ryNxDEyEt3XKn3j11+d*PA)SW(kfymCYBo7Dh0r*(
zVB9E>6qE302evB<f~J~HFy;%(1JLsV*D8deZoZ<$9NvQiR4QjV7{xPMoC(H-UZ}AS
z)!45k3nsiCw=3p_<bO)9$8hj;<w&;W$M6--@QkTq!wtCa#a1%ZaS6xIfYo4vjtl=-
zSqwNOaTEtTbi-;lEDwfdOmpLj`hu@H7hn=}a4;;T`~fpB*^pbil!1r#i+x-nO0nDs
za2_07gB^53o_@0fYbx@X%#a9J0wl|1E<Otdq4uqBmumHi>VIb>BoW!?Q9r-d=TblM
znbc4G57jSadZ~RXU6H+Dx{xo4UapBgVx?jVo`XRFJM|@N4%9K+5!5IZFS<GiAZTCV
zj%rOlctG{Uc%Eo(26nUyM*aV=OZ>Ol^(0!$IDb@`(^?Ff5BF8?66-m`o?_?8W?>Mm
zF)%`k`1C}R!GDX4iM)sxe}SR$GDCq06!>JE%*D)<e3ISYF2-2Fa;}1G?->SVxd!Uu
zGoa$kJfdm@dF4x)>LX<kmsuln$&^{}{j_G}0V$#Yeu551lfFQ46ssv@G*=f4-23tv
z+|_A)HCI`XVjZUf&nB2-*AJ9wu@rV#WLw~YApfK#v47yZd}^*x*|K=b<yaI*>1O?S
zGch2g0?O8CWh)M>W2_RRw)SN^IQ_M`C02R2&uK6aoI<j^M({;iGNmyXMsu*p6h-p)
zEu`0O<Gsyefs}{79z8KOy+`uUxE|a9$U!mDKe$uaS^%Dh8kFpj`p*71RSl=@eU%8R
z9t%qv!++Akr^aO5D!Dzk`?1FQ%of}E%#>`6DI5Zw%Hp+JEECSh6aGD}FeoG@QyBv*
zjpkE<qkyJzd<_i-)Hzl08Yj@ca(H2LP<NKA1t+d)%W`o1=$eKKmVsV;(~Z@VzKMGl
zEGS*&j&E$oj6y8S6A7k>d86I=P+{?Ap%00L4u4?;*xp+T!s4I$eYq}#Nt_dc-g`?s
zG)j7NTcXuI=YfD44v9!a&7}-VypTbK#}o=O8&kv!^L42NlMED%GM-eo!&sFEI@hyy
zDpLuyrdd7A#XTh?P*N@AfrhF?@~t&Q+!iqtYltKgXNYTMrz{kK2)K?94>P`HK_jdb
z$$zPOT&rrTyOQ7*dOcxg<?xf6I+g!m@x4<EETNB;!)ciJ&Y7o4f)4^u3uS>t8!Sy<
zR^<Dt_CU2ZIuqw(CLd|T`uT{z4sY=;>6ZRqAO8g~IL%~CM7!=ivuRFlqgFKM+BwPz
zHz*t?g)NK0VZ#A#{ye4Yrr3rK7i_DU{D1z~O$KU~f<z5gI&zZ5Y>yEVP|P6cK+oQy
z5A;|PQs$)wlD8_B5=tk828JSCBRvfFf#WHX`x!SWW>S1IrhVWB>YDRFP{*oj3d2O|
zxLNkydD1b`xu!<U&mkEP5JDUf(sU<}^V#gcv2Q|9ENYAAkt+Fm@QBEq26=>PtbY!_
zTk+;;j7ELwK*IW$iqJ%Y5Y3E8Vf!f}h5dYnF>zu?|K0ns#Plp$?9|uozCKx4)UJ2j
z`C1!7woy}J2zr9&v!Tt|*k#4#TG>QirjVjSbev!>)Z3j*&7{8AwsONpZxOA-BRwJA
z;Fc`9*EJu43DyPJA@{B)I0il<lz&B*i4<GY{MB}sL%wl?Nk>Pc4&4(^;Z#Iu$eoQt
zVF#41);^Ni_92uGq>~;mECRex&v=;%tOUk5>8vKWSr|DW<BLd{v9gOU1M)4UWIPR&
zGUs|uxU^5Wv`?0_p2ID)eGc?Y(lN2BD42qxiTv4#H8rs1j#c&D3o2qX)PF;HJl`6U
z$x$-)f9Q~OnI&J3)Zsjx!5f9~+_%&w<{JeG-yO*tcrGVD1TQF#z|dgCoXy)eQ%t=e
z5SbfrgghqrWr4^dX9~P(G>xq~CXqv#taBRj?S!sQdI#?4zMNg>gt-phR|VaBe=X&_
zP7gM7TtK@nv%~m7XfwRSIDh*a3s?C|u|)<1`L2=DLN<mMJ`v)@V$MaegMQ`(sI`QG
zCU+U<q?kqrD75j#iqS64R{h?Yb9v`p*gzbP9J{e>s7Xp-4bQ2zF*fj(C%hRu(5Uyg
z;NxNyj9UFL#$xiM3l8aoiywjGp_1`P3@68fdCW7s+!5Td<M5xUPk-kn!d_a(yu`%q
zWFsj!VCN{Ii|#A2%d?WM1n?b2QhaC41-rnS(S&SY)1j7NuUjpf*kpU%BnVK|zH1Mw
zyvv}_(Biec8$LJl!n;CGQ30=0JzW;I6hOT#VH?v-r0ZYcE~v+aPwv=>Ez3^h*@->G
zPHc5_oNh6B;vPO&?tfv!anJ4SIkd{m9O^zFU01kw?a~^b-<%qsbx*z?+H|_Yo{sXg
z;Qr=wGC#uh-eH)43R{chdg<nX<daM7Kh_W{D^GN_yt*paGH#PWlRFQkUxi#g5433+
zB5feGRREeZuKel#ed*HSwNR>j&8343E*+LtN6Hk>Uk{|p%YQ&>0c029e!pJfQUr6~
zU%-*oryV|CM(~K3-)?t2$7ufa(g4GC5v`c6-grY`-@LeOfP$~>C7FiPtn3w8z%B!w
zw6}6l*h%4ZlMv>on`WBnRKIdpb$QnU0k{yK9~QT@=!fv*l40~7!FnAGl+AoK9Ou(h
zk)4V~mCa4g;D4OtmEJ;%0e;2h7r0<UjVWybhc{NIfts@q;68)l6C?;tcjM0kNDv<W
z^7i-Cn9^3-i~Jt9)R+A_PQkTgS)Tjy*7L@feX*;Yp+V@YN1r&@k8q&ooW8)QR}*K<
zA7`8VZeRaYkGnrtWxs|3<PfUFB>bmex_D+I3eI`(secP1*day=j~;$O{*5l8WB<?{
z-3YqC&w|5d^!3a1=}Rmf+~V9ECdi*B_V~*yfXaXPWEu82)Wdq>%i**_IU&5mm{3(N
zDY)&cMUh>79IMvFsk7J0KvpC|Ieiz^^@pmDzjP=M_$^ak4Z|X`JjExd|7R#9Q1HVf
z2(tv2@P8l5T`!5+fu)fs3R~FVvhJ;e!4OX&PM1V=?hT~~L{eNJFONiZO+(+G9Wq<H
z9U?3xkFk`8U@bJe3so@};I*U^6V{E-REJQ#ik!z&UQqd9M<{$@$9OgguR!_t<vzxp
zi?zDWo$bEfWB<LEM-OuMYIAVNl3&{`mnTix-%nM`{NJ-7gP7qe7q_i;MZUc+@)*z!
zDsqpC@3gFQna=xi@SbAU+>07`z0_S)-WA=QboA0faY6@W3@9#yf1w!nDV&Lhe06j6
zA1EB;lCw!<B|CrHFt}(xMILlA$2kGusgU{P94DuPmn7j#pyt15$j9C^;G}1=n<eW_
zJV;ni$(rInc8=3q%~M?I$;7oI$0Rg9)y0+r7hBGBv1PuCE&kXm03HL~CrLbei|T>s
z0XlUO>0R#8qx(m2l}76WRPu+QvU^ck?33P98ooe5HW`0tGL|pn%fJ`55x_cAKKxbX
zPdA&}yXF1;fXc9Z^M3qpz5elj{ORiYhnv;%vv<p@yVcF=?Pq5C#<G<8<SNJFcNpT4
zcb3OweF}q5ZkDV~xTgG3JxNErlLYQbd-U{kHt`O8KYuNl!#rd54?X<7VlHTVvXA@!
zqrxGj!GwQ#vfh}%<(aE+r1}x=92K@4InEr}6r9Vm{}gd(5@Sv!9Z!9GWz505XJg{I
z>@o9SKkvytp5`4D4jQ<n_$f}p=S7NLq11T$_WRqrx8JFYhxLqKF8_0Ly?pi63;F%)
zBE~o*s2K*N5-cOP(&ao9{xsR{t6cH$?uZ?Bu-kw3z9-u+7Z<r$^mVQbPDveYOuLe5
zWn%xsl0O8vR`%`ANjrZF6YI^}vLyDBGP%fgNExg=tz|c6U#6fZ|3+hLetSv9V4bz>
zdqqdB9F|#gZHJ|?ve0(PqN`To$jbkw|3}L#Sj;J#d>i{=4w<o*S#)WP+%&a|Vq235
z_40q%ot6%j3w4M(>Z#Ir5p7zSPWQ;A9O<P&$Tr+z%Y;ik2fDJ?EBlhJAZ0-|n*|(A
zRW+3g<+COHRfP`5uL`0K2-mg&Z&`wty;3|`$6;b72e1;aA}jzi;wf1{$yiCDi<4cF
zzA<?`PUyhz2-wn>aC*8H$TtRl1VIA}R&#$_#&HrP#h_~7(-66n0F@N^4K~V;rp#9;
zA2l4Aay|=^jw?)t{PtWW)>rt!fml)Ikg{cO0@t1S4T3NXt>9H)zNedQA{V@9R-$Da
z{KUZ}B9%0jiLv;DoV4+y2I<FNe^kWTnI9PNX@n2qClXPfvL##Dj^PIymeRpfC4ztZ
zw5UvJe8K&P{P%s+!w9<?U+k3TEi!$ep}sOWg;r3<%DY#1pd5t@qBVFxz!AvJ?2Bkn
z({d6?mH`CDseAGab9|P9$N^L<m_Q@hKKQXB-iaE4R%I+cjm=y4<xVTr=Mqg|Vhc&-
zWJJ)1nkJ2X=1P4cMv)%3+7n1$<&b~Rm&}<{7LH&H;Y-lmv-gp>2AYKj`4BxK%_Stk
zQX~@?RD=O&5zcUCawh$X=KHAhcq!Z60y@9jj6tP1-nyM{kBj?ay$?PrYEx*ECp<sX
zp%w$BtP---8jriEJc-REi>*e%GLmy*ht(e@BUWDMEOPa$QdS=YcJgZ;`bvK*t|1F)
zQel{lMFO4BD|;Gxl%U4}q8K%oK&q8orxND?f*|}_|EpoISD6>$AC_C~_tc*{fZBg&
z;u4ZbyIx%f^%nFyb@u_Sj^yfYk2n*$i58@G@I0)lXe-o!I*RJqUrWTVHVuVqt(WSz
z8an#Ovu=6kdcWQh`-s}qLl1wbD}>S-xr^%8V&KM7wRgxhFVO9iU)mKRjEqy;q;DI&
z=5M0#HvU6?W*?>xTU2)n+im$E2*_Z4%24YadeS<<)b7T}ee(GrT_v0(h_1+Z6_D;}
zV^WuV>^RZ84v-P2KDbI$ENqHjaWcf!dJVu&DeE}3%r;5@*_Axd<W_%AB}kBC_)o|j
zd`nij6TmC^3iT11q?6zid|-fsy(3uhS78J=<|kFAjVvRydN%AOL7>jk=9SX$>&Lz!
zO62U7(dv~YA~R7{jV?6{?CMTjPA}7xUE8BZ*_4#Im6q9b@>Rwo^1~sSBBk5TA}Gr1
zxXOwNpC0ulJJvd6clm$+YW#=s?{{weyA8CCLI0kOpC+DZF3!^+YP~AyDP<j}mZ=mm
zsa@$av-j4Lw472TZ}mLdq>R8V&!QX<vaUG-<y@-g03*Gb<66l=?RK2nub(l9Q$fEr
zdX<wB2S+)icTG_S@72)i6-<$tC<z1mn#xFxM+s+q_q3^(ZJvLJn|hxgmbH_t&%L(t
zRiWp1-b3o1)D6_SlmM^GNpRP_;((2Zy?T_hWDZ3Rgme_yP({Ug@uOLQki$W)?UiOx
zUNrDnSqM^gL;b*+B3kmaNVX9}PRK>)C}yPPRF=FLZ+<vd8lAk1&vR~+6r7K-_h=;%
z-<8J5`;2z+br3TndeNc>!~BowVxmrfQxcm^Ewig1=+=}W$(uDNp1{?72OHKrUj%U+
zfmU<5_#c*ageR9l^#T)vC#HudrU8d1rUJJorUO+@1T-};G?#$%10Mu5H83={qx1uW
z4N^8XIWsv#G&46tK`}u^G(|=+I5sddK{Q1(MMFe6LpMGkJT^8tGdV;wGdDv)F+oK%
zMMg0=HZU_mG(|H-Lqs`4H$GhoFHB`_XLM*FG&ndQlMz`ce<jknPgGG9$MN%fcP=n2
zGYBf^Acz|ZE})En3kvQ!io4^23WU;zn#!11-Nx2f(^y#OA3$PaXe{Y%EbKHm2^I1^
zHlO_NoO|DUXU;V<+sMq~eRu#5!DR1;jU>n*FSDZm79h7t+5w$#0D53Qbiq#cMW};%
zXrNGnJ<!K*e_NnYmY?y*H&NI}YK9glKoMH0?j~)AHd*_IA8UuLRJ);r<u=kT=w<02
z)FJC$_7i+~^cMD#{GlaT|0l1~LgZ}Zd6C{CUB+HndC|AfLiDr#xpW@;p-(pc$G04W
zL%hS&F&L0l{{Op&;RqasVHlKE*S$If$KeE=g0gJ>e~nj9!U&AQ7@UW3I4xTk@_lFE
zEL?(fvc*?kn1D&R02kpFT$a5Z@hw;2DqMr>a04oEleXwB-h(OGyVri~Hr#=`@Q8hk
zB34O0sK@YxB_oSZ;TcPx>NGrO*+zN+FJT5|VUDUP)T+EPaZghwTblA^kbN(B%1L82
zPf=Q$fA>@;efjLEURwR;sX_X=<|(giq3o$q*~br_nv~Ui=BZg>+OJP)QTST&R8Uy$
z@>Ep#QSsEO@N2~r*E?}7(f}~?#LN>jPuz$!2#h;1?!>qgmm!sbNhc<qMiu^cnsV81
v`v=eWUg8R6Ze(+Ga%Ev{3T19&Z(?c+b97;Hmz)Fy2?{tlHwq;sMNdWwY+AxK

delta 5225
zcmV-v6qf7Rix97g5U>qd1Tr`?G?NioDt}vZ<Tw?6&#y3Vonl&J-7E`g9$<Hef}(bT
zO%+f(`{3zT60_ZQW4kAlec-oqj*jHE-0mbb1d5kLmUTJz%jeFsyDZOsca?|#zrDG7
z`<|Q3GFyluySdFc%X1Fz@{DI(N@j{OE3{!k>g;Bp{jl1*s<~cEDOaDab+M{;M}J3;
zbx;3)c3q1jv8tLaO%3;Mx}~eDW;5H>HjXxZ$aU^MT??~v*L>AI((6sTJMP21b(5Sb
zR$JHCcg@>dR}II`UH^FVhb(8doO7Wx6J=nc8oK&3T*cY-S_!-Q`*o>T?eTgoiemL}
zElPOU!TnpDG5#vPYUrV92U>UM!hgV1f`h)=hf(MH<8J8Rz*EUrecf!}07nPF`F0i7
zZR=0<HjF>wy^pgaL1ZS1j4kn|cviqaDT0D-{LV6mqpQ2Nsdn_VY4O%n0XR>hpt0Mz
zJuwP|3_VRg)Vpv5(A3cABmP^36F_yyrWLD)`+DcbwKQnjW<64Z<eyoW;D1aSq8*M3
z_Ru&NjmKM)@^_+9_y-bOSP|p~hrfL<<kXl&&Uk?r-fZD?K}#f`EKwj#o|Ukyq(ui@
zboI^6)jzKyrJyq=xCBZqGeGKQfAzzUdA5a@s0|6~c<@X2;kZ28W#3=@ncOR(=yVh@
zJs~S`!;;BdM7nR=E;z}tBY*aLP||^Pu&F=fQf-2pR6Q99Tq`EWt-0dND3h&)W~D9Y
z+^^dG;W(f~%jiv|KpQ^n!Xy|C_}8cFLM0}+5HoC$mCOsKfpl#(Q$~?mzim-N4@hH<
z)=z2#bF4Z?lXq?1+|khM8{QrJ@Gkf{m@3*PMIW$d*Q^)qNtcpoZhz_K!FAC}w$)Ij
z7PT4qpT?+>e;;;`*z>5>-X@!rKb$RCpayXl62xwqRKdbN<arU6+N@yK$P$+t#;pcT
zF$Mqp!e-$>uvD8(*1YpF06RB$mXVgG`GOU5cn<+k>w*>F6i-;Okc`W6X2#{zjN@6#
zFd6N*S#U2TS4BGp#((n#mi#EL;R}-Cf@za+3*mcqm73bP0`>)n8eA|);cp9vfut0!
z;^2p2Ssj)YlVyUH+(x3_@ulDbLV^hnj-_>EP(;m6rM1o(MCf?f=LMoPE3Aaz!67vG
zMK=`bvo&~=DN-`0M8FZCS*8m4DilnuZ~3~s)L&?RRzVRVUVq2@;#sd{e)5&fPyP?h
zFK6Z4`gFdad&x?po>RSAQhg#y%``ku7EQ!iUh?LNIZkf`Gs@+&p$-xV`q%VEy`&#x
zP@g%^ndJ)bqxEFf|Btxje_LE<*;=LbqZYwxF=am7yDk*grwV(_od;iqLAJ)k2qoe<
zN29|_!DM05!+)O<sJbjr5CSDWl@@cUFclwl_m_(~R<eTYNwxQcfJ9LObtw;MI42oV
zH-fT?t&ICj6~txH$U-q~XL3I*Ir)SXX+WM}!^x(1s1CDOLMD4l!yvrR+2Em1>)ff!
zf)Z;#Hh8|k9NMm@N{i)i!>rnZ2n6*f9f=L!Wv;nKV}BFrmdm+lqNST<_D;lrmJGD5
z$;n<F*vEJyMr-Y=rVsY(YfJ3%ZlAJXAUKC&d5PeQv<O{dGEDYhlPS&Wuh&pshl@A9
zjs;R4x_XGj`0^gfL+5%33!nzYLVp{$@Vx*$_ca*VBh4KnJFW(9yVxaysi(?P*0S8h
z(wLlEuYYcj&1R^vKl9CYaT1DH5(<Yzr$u^JgKffbf53l_3j&%FBeY=PrOAG3a5T_V
z4KJa=f;s0pJ>vk@SM|?a4(84ZJ(I*6y0YqnJlfrk8kU}3eA^CgPTwd!4;ECRbH_J+
z;)FuX+Y<$@h;^gCdF!xw^U#OFMu#v09B(`Y(|_il=3~1qrCo&}1he;^b{JIj=C(qw
zeJTP0H5>|&gqlkglzgUwnzAW1R5sS62lHjEOg33)8dbWfZu-GhTjJ|UKUKMcUehiz
zbMa0E1(eb=eW0Uix%$pKB3_G_iFZU&NpQrawo@5RAOfBf)58hh@}LoRiWF2MuH|;_
zu74u9g<g-CS=Ij(wvO#TczhSk0$b?43pfpn(OE>A6vQBiw9poKw87K#wWHYQnl1I(
z7)(Np8DpeR>&HF*4R}ioNw@U><@sNNfODZzBAVssDaJj8&sxczHz81txPb|n6plOw
zhl2pPMSiL<O}PyNF8J1NjQc~otVFb4U4JG4Ql5NKrzpS(fc60dAa6z|-Ue-w5Ce*U
zieEx!@)VZ&I-q6s&yHX<!qV&sn!)ioG=uZ}gy-SLmj1UNhYIVI<gi0mH=Fw4@j-vy
zejS@m$X!;C6%q&xfjJ-JgRw2-zp<n&wJ&KE^v5bUCjiKJ0LbG&gdLs$Ab9S;)qfvH
z-B5;Wcb`U`m%$BA3B8<h;5FRt2@~_gde8U;Cf*Y22S~i()NO#5`u5Cgs2El`RGE<A
zs2C0eC<Qw~u~XIHoTmi`Q8h}L6KXV34bBJEES7nr2PjUVa7x7m#Ax(1aQ9L{#RN&s
zYV$Hr)6hZ}MS2f}KV$7!ryyUxPJh9{v8Uvea2-pBCvtUkEm=Z!ZX<CAwyQyocg)py
zPZP>9*F$hm_Inf2L79*u`xZ#opw-_G^dbZZjd5W(b`8yv<-jD|>Y9_1eAVl{7#2|k
z)@582!!6ap{>_+;E+`ShS1D#ZKchrtIOLB4_LL9;#fy}vOu=+HlMP#pY=6=irop=X
zlJCY4Hu@q2tNSW28)iJb*@u@8lj%+T0u&(V4qS6TLs3H!g`t{&U3iBn?N=Tdir=4R
z#Rim?Rw<8;SZ-p?!Od*XMN@&&yD-q-Y3LwCsJ+Lc+m;i+*G6N2K}d_;hD>6p?!pV7
z^QJ3@cJRT*mQ(KzMH@H{4u7iyeuUc6s*iZ!<Lnen5?w!}YKo*QXnTrS(N5x<Peb$g
z_LMBwge<H%{1?TS%Tj>-D4EI<Yox;tD;fwLAsWCGeHz56C`v-#@f{%u_|E$mPC*1f
zpYwCO7Oezl!)am^lb;QnAO*TjOI_x<><|nq%NNS7f3fiN%Y--mEPod$f_&`Z|M*?5
zp%fi`3CCE`AzlA~_(C%wLGvI^9C>ltC{COqapI^YcG}tIQF{2`rH6wcy>N4s(2_8P
z0nYLu%mwMa2~}z=SIHbOQjA>;jEOtrxr;T?_^ISIEjNMD-eYQkX10{8H}fzDBp*T*
z{kevO>>|-_{*|a&Dt~URfFVyK)R%*V2&bzQvXhHgFQf(?x%<~8eE-W(%NOAL^8!C-
znisEvSsEbQOz->Uf|L?O`ThcAtG|r-^ccw#Vt#$S<^|RW=lHyZHo+t(=8GSfRgrwH
ztDgg$m+n$fOW38m5)W7luu1;`40<Rf!?#8PbI&goP%-A^e}8_fL9A=aU_bjXig5^k
z$I%ZlBKW9-iE?-P3QJv56SBtG&*pO0(2Bfd+=N2sD=0A_u7vsm58TmUEL$Mp4erp>
z@Du}htYG+vzQHhjT;~(?4IcjV)o*Apmb!Eo#Xat+uexPOfZKwBR$@09W%E!&h&&6G
zCln6bj3E{d{(tu_o858zu|vIBY!_gVC!hR&Q~y*C>p!}xTdE)nNL6xF{`rqeg0Oi>
zEyC$nB?uNz;YfIlNC(Pq2pJvc`*t5z(6@3iVOG!&O~zlDq}stNj=SxMu7iw`{W;q1
z;%|(+z#jX$UoLz(!sa|E;XURA=R#3%+c`J=MF2wggClD~o;pUWoaicCGG?JnN{|)A
zxf1e%3c1QAOTM|e`Zv6k@slB&6SIV66FYwcOun!A1bNWO9Ondpr$XkFbDW$GUXp|}
zftvrKAs>6wfRmodZkDV!@gQM6C2NZR*f~ybHBWJ+Cll9>9Fx%aR2N$gTx>bh#g_Rl
zw)kVO0C)^^pCs|@Evg5i2k6vEq<6VTkM1AART`}iP{|*H%I-yFu}^waY4`#K*<^p9
z$ymOOF9ToLMgZ$f`S4ejKizC@@0R!X11iJv&HM4Y_4>#A@u#cnA8uC5&)zMs?p8Oe
zx1X8m8_QDWldBwy-(iSH-dP@#^(hQKxmmI{;hOSC^&}nfP7=5$?a|ZE*~B~W{rt6H
z4)cuJKlJeTin*Ze$v*D?j|zvB1`~hg$$Dc3muIfRk?Kddb5z)J<T!I=Q*bWN{!_%E
zNsKv_bUgL#l`#kJo{fp;vd7GS{k$jpc$#-qIB4LO;-@$XpBE{1g;L|~+wX7h-hQVp
z9@aB{x%|(~_43tMFXZ>Hix}gOpk^45O0bOFN|*Cc_|s&&uX4r1yCZhk!ES%s`<`sS
zTwLU0(bu^$I3;zoG3`pKm5KciOa2hxTG_WdC+++#OsqF=%aYhj%H$%~A!V@gw3gkN
zeVKxu{2Ps}`RyeYgLT%j?-d=ja#&`~wH=nm%0k;Ei>_LUBP;)#{vR!~U@@m`@@?#g
zIb_CKX3?cFa?{i<ifv6M)XRTkcUn4BF4Q6FsHaNfMYL&QI^83ea-^3AA=_|=EfX&F
z9O%kkuk1^@f|Lc>Y!+}dRn=4~l+Tv%R~0%Kzbc3}AY9u9yk!Ym_Db<&9fygT9KcGv
zim(98h^J%)C1WLpE>3ny`o`q(IH3c-BVbEk!s+Q+Am14H5d;k^Sj~TJ8OKSG6oaaP
zPebHR0#s7uH`pjYnlfLZeAIAc%K0ovI<7Do^4oKjSYP1>2VzB;L&}!D30!yPHweNo
zw1QWG`JQgJiCplaS&5cy@Dm4|iurAa#xgM$e~^<le$*iS`0I~~HgN91rx8AcpGZV`
z%9d<pJBA-<SV{*^l?Z?C)1oq^@dfuE^56GO4<qbqe6drSx5)H?hWg6j6k0(YEAL+6
zfpQcsh}Pf%0Y@M=voE4SP0L9nSq2aor|!u!%<)+YA_q{dU;>R~`{2ilcqeKET9vW*
zG&XPDmpiRgpG!1>i7h0RlMz85YMM0mnJe{)7)5&AYEK}2l|z3%UovM-SvZ0*gfBsJ
z&)!Gk8fX?C<U{m`G?$PBOOZ@qP!R^8ML5Hm$(i&kn(w31<E3nO3+ViAGX|C7c<Xk)
zJudEx^*;Ehs7;|sp78ukhguAjvP#HSYdr3v@+3BwEVdd2%Sg_N9aevsj97W0v&hx2
zN?Cmr*vYSX=qrD%xP~mGNrhoH7727huk2~)QGy-^h+@=S0;yJVol2Yo2!ily{jY|-
zUS(d0e^_p{-&23;0BZl8iAzW#?Rs?`)LYQ+)ZGWPI+ClqJ>pF0CR&i%!Sk@HqODK^
z>L{vbe=QNe+B6idwO*><YUt=A&${KE>-~C5>?3MZ4?TaNt`JIV<SwdTi-8+U)!rf3
zyg;{4erZ>PFfvYUlfG^An!ky{+xQRpnSGc(Y*F1QY`5iuARvSFDMPJy=t=7YQ@a}@
z_sQpjbd_+DAi5&sRY1C@jY(bdvExMVIzUF8`rs;2v9KwA#mNv?>oovBrL5!BGTSHt
zWLNS)lUsj1l^{Wm;XfgB@GV*8P5`guE7V75l1_q8@PPpi_Ksl1Uxg9in4eUcHnNP+
z>e;ZD1c5qBn^#K1uOIt{D3P;QMypqrh|EM)HM-O+u&X<9IlW9%c5ROuWm8h-R$6A$
z$yXVV$Pb5Pij;0Ui=Zg4<0>mAe0tQI>{#oN-Q|D(tMMPkzu&p>?>5jn2K{?Bewuiu
zxj0XQsP(F(r<8S^TBcIOq;{px%-&l|(sD|Xyw&q)lQIIgJd1Kb$hzhTlyj+`1B~=$
zj%y_gwcBxOzkbFbP6hqi=v7Wi93174-Ze!TyjMf3S1?6pq9hFLYbql(9wnUd-P5LC
zwt0UZZt8u0Sk_LmKKI(nSB0M6c@L?3Qa4cRQUbg#C&69!iUT$t_Ucj2k~tJP5Yka(
zLlqV0#gAqILJkMHwpW@(dC|aQWg$q}4fO+OifGBxBH2a^IUyIFqnMGFQ(5w2y!qi+
zX>{^3KF_&PQgA-T-lLU7d{-JH?=#xP*FiIk=tYYj4D&yxi-|e`PDyMwwal)5pj%Uh
zByZN7cmh}R9c)<hd=bQP1X|7I;(zWOgerqQr-wbK0f#-O0=GS<187bJGBYzbm#Osw
z9|STpGdH))^#g<rQ87V8LqbJCFfv6oGBPnnMM5z_Ml~`<LNP`&ML0r2I6fdeF+oH_
zLPbF^GDS5qGBHL)LNP%`H8Ms*F-9{*I6^`=K3xhgOl59obZ8(nI5(G}0|Y659n#NJ
z%~2S~@#phB_td@J+r24-6p2!y(B&femETb$N}=K(XvU7UjX4WDV_~e?*jV@nSec<6
zdz%(kW;89V;{9x1^Euz=obUJh%yZ1lHZrq#2d3aInCwl_NP-LsGOPHr0J;C9?XV5H
zVHfO!ov?#_5$d2G@)X)&5A^VVY!fud3Jd=FB89!A5;Q^?D$qoAH)$(0%UU-4SPN{X
z+5xRBw~*SQlcoPqtE}UbpWw%%x6noMm$u3JK6sTDB4;D-i}V)hGIq+U&wUFmL_h1l
zOXsl<dSoMie9M05<r|&`pkFpN<kbUk5Dvi*49ceedi5|IfunE|s<J14vtB(0$Dsx%
z;4F;5DcMrq_YK2oxCm!t&qlp)4$i{`7=@c~N%o@bTQ0*DxC+<cI*h>$+M>5O3FER?
zGk$CWZozH1$38|8r%67j8Mx1qk;Pe<W9d_U01sI<lODl5EWjc>rfLf71AH=ZPg5pa
z?)GJneOvXElWHrTqO?3{?5R%r_|j9o^y96kytMk&Q$g8M$y0-}_fI_)m90;DDk;pp
z_SC5G*@v9U3g6c}RTO@9d1_Ml_0AL5J8>;iKQQye%o8(D+=w&?j5{&z#JCfeAyt7%
jCnlY03TwYjIeoEz-y&Y_m#_r{9||@wFbX9lMNdWwi%|{{

diff --git a/doc/src/DimRed/DimRed.do.txt b/doc/src/DimRed/DimRed.do.txt
index 0ad06f689..0b7bc1eb0 100644
--- a/doc/src/DimRed/DimRed.do.txt
+++ b/doc/src/DimRed/DimRed.do.txt
@@ -782,7 +782,7 @@ Note that the function _multivariate_ returns also the covariance discussed abov
 import numpy as np
 import pandas as pd
 from IPython.display import display
-n = 100
+n = 10000
 mean = (-1, 2)
 cov = [[4, 2], [2, 2]]
 X = np.random.multivariate_normal(mean, cov, n)
@@ -858,14 +858,14 @@ plt.show()
 !ec
 
 Depending on the number of points $n$, we will get results that are close to the covariance values defined above.
-
+The plot shows how the data are clustered around a line with slope close to one. Is this expected?
 
 === Diagonalize the sample covariance matrix to obtain the principal components ===
 
 Now we are ready to solve for the principal components! To do so we
-diagonalize the sample covariance matrix $\Sigma_n$. We can use the
+diagonalize the sample covariance matrix $\Sigma$. We can use the
 function _np.linalg.eig_ to do so. It will return the eigenvalues and
-eigenvectors of $\Sigma_n$. Once we have these we can perform the 
+eigenvectors of $\Sigma$. Once we have these we can perform the 
 following tasks:
 
 * We compute the percentage of the total variance captured by the first principal component
@@ -883,27 +883,35 @@ where $v_0$ is the first principal component.
 Collecting all these steps we can write our own PCA function and
 compare this with the functionality included in _Scikit-Learn_.  
 
-The code here outlines some of the elements we could include in the analysis. Feel free to extend upon this.
+The code here outlines some of the elements we could include in the
+analysis. Feel free to extend upon this in order to address the above
+questions.
+
 !bc pycod
-#Now we do an SVD
-U, s, V = np.linalg.svd(X_centered)
-c1 = V.T[:, 0]
-c2 = V.T[:, 1]
-W2 = V.T[:, :2]
-X2D = X_centered.dot(W2)
-print(X2D)
+# diagonalize and obtain eigenvalues, not necessarily sorted
+EigValues, EigVectors = np.linalg.eig(Cov)
+# sort eigenvectors and eigenvalues
+#permute = EigValues.argsort()
+#EigValues = EigValues[permute]
+#EigVectors = EigVectors[:,permute]
+print("Eigenvalues of Covariance matrix")
+for i in range(2):
+    print(EigValues[i])
+FirstEigvector = EigVectors[:,0]
+SecondEigvector = EigVectors[:,1]
+print("First eigenvector")
+print(FirstEigvector)
+print("Second eigenvector")
+print(SecondEigvector)
 #thereafter we do a PCA with Scikit-learn
 from sklearn.decomposition import PCA
 pca = PCA(n_components = 2)
 X2Dsl = pca.fit_transform(X)
-print("Check that we get the same")
-print(X2D-X2Dsl)
-
+print("Eigenvector of largest eigenvalue")
 print(pca.components_.T[:, 0])
 
-
 !ec
-
+This code does not contain all the above elements, but it shows how we can use _Scikit-Learn_ to extract the eigenvector which corresponds to the largest eigenvalue. Try to, based on the above, to address the questions above. 
 
 !split
 ===== Classical PCA Theorem   =====
diff --git a/doc/src/DimRed/Programs/PCAsimple.py b/doc/src/DimRed/Programs/PCAsimple.py
index 9e7716348..c23d51e47 100644
--- a/doc/src/DimRed/Programs/PCAsimple.py
+++ b/doc/src/DimRed/Programs/PCAsimple.py
@@ -24,13 +24,9 @@ print(np.cov(X_centered.T))
 x = X_centered[:,[0]]
 y = X_centered[:,[1]]
 Cov = np.zeros((2,2))
-cov_xy = np.sum(x.T@y)/(n-1.0)
-cov_xx = np.sum(x.T@x)/(n-1.0)
-cov_yy = np.sum(y.T@y)/(n-1.0)
-
-Cov[0,0]= cov_xx
-Cov[1,1]= cov_yy
-Cov[0,1]= cov_xy
+Cov[0,1] = np.sum(x.T@y)/(n-1.0)
+Cov[0,0] = np.sum(x.T@x)/(n-1.0)
+Cov[1,1] = np.sum(y.T@y)/(n-1.0)
 Cov[1,0]= Cov[0,1]
 print("Centered covariance using own code")
 print(Cov)
@@ -39,24 +35,27 @@ plt.plot(x, y, 'x')
 plt.axis('equal')
 plt.show()
 
-
-
-
-"""
-#Now we do an SVD
-U, s, V = np.linalg.svd(X_centered)
-c1 = V.T[:, 0]
-c2 = V.T[:, 1]
-W2 = V.T[:, :2]
-X2D = X_centered.dot(W2)
+# diagonalize and obtain eigenvalues, not necessarily sorted
+EigValues, EigVectors = np.linalg.eig(Cov)
+# sort eigenvectors and eigenvalues
+#permute = EigValues.argsort()
+#EigValues = EigValues[permute]
+#EigVectors = EigVectors[:,permute]
+print("Eigenvalues of Covariance matrix")
+for i in range(2):
+    print(EigValues[i])
+FirstEigvector = EigVectors[:,0]
+SecondEigvector = EigVectors[:,1]
+print("First eigenvector")
+print(FirstEigvector)
+print("Second eigenvector")
+print(SecondEigvector)
 #thereafter we do a PCA with Scikit-learn
 from sklearn.decomposition import PCA
 pca = PCA(n_components = 2)
 X2Dsl = pca.fit_transform(X)
-print("Check that we get the same")
-print(X2D-X2Dsl)
-
+print("Eigenvector of largest eigenvalue")
 print(pca.components_.T[:, 0])
-"""
+