further update, added also

doc/src/DimRed/DimRed.do.txt

@@ -703,6 +703,44 @@ and since $\bm{C}[\bm{y}]$ is diagonal we have for a given eigenvalue $i$ of the
 In the derivation of the PCA theorem we will assume that the eigenvalues are ordered in descending order, that is
 $\lambda_0 > \lambda_1 > \dots > \lambda_{p-1}$. 
 
+
+The eigenvalues tell us then how much we need to stretch the
+corresponding eigenvectors. Dimensions with large eigenvalues have
+thus large variations (large variance) and define therefore useful
+dimensions. The data points are more spread out in the direction of
+these eigenvectors.  Smaller eigenvalues mean on the other hand that
+the corresponding eigenvectors are shrunk accordingly and the data
+points are tightly bunched together and there is not much variation in
+these specific directions. Hopefully then we could leave it out
+dimensions where the eigenvalues are very small. If $p$ is very large,
+we could then aim at reducing $p$ to $l << p$ and handle only $l$
+features/predictors.
+
+!split
+===== The Algorithm before the Theorem =====
+
+Here's how we would proceed in setting up the algorithm for the PCA, see also discussion below here. 
+* Set up the datapoints for the design/feature matrix $\bm{X}$ with $\bm{X}\in {\mathbb{R}}^{n\times p}$, with the predictors/features $p$  referring to the column numbers and the entries $n$ being the row elements.
+!bt
+\[
+\bm{X}=\begin{bmatrix}
+x_{0,0} & x_{0,1} & x_{0,2}& \dots & \dots x_{0,p-1}\\
+x_{1,0} & x_{1,1} & x_{1,2}& \dots & \dots x_{1,p-1}\\
+x_{2,0} & x_{2,1} & x_{2,2}& \dots & \dots x_{2,p-1}\\
+\dots & \dots & \dots & \dots \dots & \dots \\
+x_{n-2,0} & x_{n-2,1} & x_{n-2,2}& \dots & \dots x_{n-2,p-1}\\
+x_{n-1,0} & x_{n-1,1} & x_{n-1,2}& \dots & \dots x_{n-1,p-1}\\
+\end{bmatrix},
+\]
+!et
+* Center the data by subtracting the mean value for each column. This leads to a new matrix $\bm{X}\rightarrow \overline{\bm{X}}$.
+* Compute then the covariance/correlation matrix $\mathbb{E}[\overline{\bm{X}}\overline{\bm{X}}^T].
+* Find the eigenpairs of $\bm{C}$ with eigenvalues $[\lambda_0,\lambda_1,\dots,\lambda_{p-1}]$ and eigenvectors $[\bm{s}_0,\bm{s}_1,\dots,\bm{s}_{p-1}]$.
+* Order the eigenvalue (and the eigenvectors accordingly) in order of decreasing eigenvalues.
+* Keep only those $l$ eigenvalues larger than a selected threshold value, discarding thus $p-l$ features since we expect small variations in the data here.
+
+After this we ask ourselves how do we prove the link between the maximum variance and the feature reduction.
+
 !split
 ===== Classical PCA Theorem   =====
 

doc/src/DimRed/mlpfranke.py

@@ -0,0 +1,68 @@
+# Common imports
+import numpy as np
+from sklearn.neural_network import MLPRegressor
+from sklearn.metrics import accuracy_score
+import seaborn as sns
+import matplotlib.pyplot as plt
+
+def FrankeFunction(x,y):
+	term1 = 0.75*np.exp(-(0.25*(9*x-2)**2) - 0.25*((9*y-2)**2))
+	term2 = 0.75*np.exp(-((9*x+1)**2)/49.0 - 0.1*(9*y+1))
+	term3 = 0.5*np.exp(-(9*x-7)**2/4.0 - 0.25*((9*y-3)**2))
+	term4 = -0.2*np.exp(-(9*x-4)**2 - (9*y-7)**2)
+	return term1 + term2 + term3 + term4
+
+
+def create_X(x, y, n ):
+	if len(x.shape) > 1:
+		x = np.ravel(x)
+		y = np.ravel(y)
+
+	N = len(x)
+	l = int((n+1)*(n+2)/2)		# Number of elements in beta
+	X = np.ones((N,l))
+
+	for i in range(1,n+1):
+		q = int((i)*(i+1)/2)
+		for k in range(i+1):
+			X[:,q+k] = (x**(i-k))*(y**k)
+
+	return X
+
+
+# Making meshgrid of datapoints and compute Franke's function
+n = 4
+N = 100
+x = np.sort(np.random.uniform(0, 1, N))
+y = np.sort(np.random.uniform(0, 1, N))
+z = FrankeFunction(x, y)
+X = create_X(x, y, n=n)    
+
+# only training data, no advanced splitting
+X_train = X
+Y_train = z
+# only one simple layer with 100 neurons
+n_hidden_neurons = 100
+epochs = 100
+# store models for later use
+eta_vals = np.logspace(-5, 1, 7)
+lmbd_vals = np.logspace(-5, 1, 7)
+# store the models for later use
+DNN_scikit = np.zeros((len(eta_vals), len(lmbd_vals)), dtype=object)
+train_accuracy = np.zeros((len(eta_vals), len(lmbd_vals)))
+sns.set()
+for i, eta in enumerate(eta_vals):
+    for j, lmbd in enumerate(lmbd_vals):
+        dnn = MLPRegressor(hidden_layer_sizes=(n_hidden_neurons), activation='logistic',
+                            alpha=lmbd, learning_rate_init=eta, max_iter=epochs)
+        dnn.fit(X_train, Y_train)
+        DNN_scikit[i][j] = dnn
+        train_accuracy[i][j] = dnn.score(X_train, Y_train)
+
+fig, ax = plt.subplots(figsize = (10, 10))
+sns.heatmap(train_accuracy, annot=True, ax=ax, cmap="viridis")
+ax.set_title("Training Accuracy")
+ax.set_ylabel("$\eta$")
+ax.set_xlabel("$\lambda$")
+plt.show()
+