version 2

Work2/report/Work2_report.tex

@@ -36,9 +36,17 @@
 \usepackage{amssymb}
 \usepackage{amsfonts}
 \usepackage{amsmath}
+\usepackage{commath}
 
 \usepackage{float}
 
+\hypersetup{
+	colorlinks=true,
+	linkcolor=red,   % Internal links, those generated by cross-referenced elements
+	filecolor=blue,  % Links to local files
+	urlcolor=blue    % Links to web sites
+}
+
 \begin{document}
 
 \setlist[itemize]{topsep=0pt, partopsep=0pt, itemsep=3pt, parsep=3pt}
@@ -74,7 +82,7 @@
 \subsection{Κλήση μεθόδων επιλογής βήματος $\gamma_k$}
 \label{subsec:polymorphic-calls}
 Δεδομένου ότι οι μέθοδοι θα πρέπει να καλεστούν και εκτελεστούν με παραπάνω από μία τεχνική επιλογής βήματος $\gamma_k$, δημιουργήσαμε εσωτερικά της κάθε μεθόδου ένα κοινό interface για τις μεθόδους επιλογής βήματος.
-Αυτό έχει τη μορφή: \textit{\textbf{gamma\_<method>(f, dk, xk)}}, όπου το \textbf{f} είναι η αντικειμενική συνάρτηση, \textbf{dk} η τιμή της συνάρτησης κλίσης στο xk και \textbf{xk} το σημείο ενδιαφέροντος.
+Αυτό έχει τη μορφή: \textit{\textbf{gamma\_<method>(f, grad\_f, dk, xk)}}, όπου το \textbf{f} είναι η αντικειμενική συνάρτηση, \textbf{grad\_f} η συνάρτηση κλίσης της, \textbf{dk} η τιμή της συνάρτησης κλίσης στο xk και \textbf{xk} το σημείο ενδιαφέροντος.
 Για την κάθε μία από αυτές δημιουργήσαμε ξεχωριστή συνάρτηση που υλοποιεί το παραπάνω interface.
 Μία για σταθερό βήμα, μία για επιλογή βήματος που ελαχιστοποιεί την $f(x_k + \gamma_k d_k)$ και μία με τη μέθοδο Armijo.
 Για την επιλογή και κλήση των μεθόδων επιλογής βήματος εισαγάγαμε μία νέα παράμετρο string που χρησιμοποιείται ως enumerator και με βάση αυτή γίνεται η τελική επιλογή.
@@ -93,7 +101,7 @@
 Για το λόγο αυτό, ενώ η συνάρτηση δίνεται ως symbolic expression, μέσω αυτής υπολογίζονται αυτόματα η κλίση, ο Εσσιανός αλλά και οι “κανονικές” συναρτήσεις MATLAB που τις υλοποιούν.
 Έτσι έχουμε την ακριβή αναπαράσταση της κλίσης και του Εσσιανού ως συναρτήσεις χωρίς να πληρώνουμε το κόστος της subs().
 
-\section{Απεικόνιση της συνάρτησης}
+\section{Απεικόνιση της συνάρτησης - Θέμα 1}
 Η συνάρτηση με την οποία ασχολούμαστε στην παρούσα εργασία είναι η:
 \boldmath
 \begin{equation}
@@ -102,11 +110,11 @@
 \label{eq:ObjectiveFunction}
 \unboldmath
 Στο παρακάτω σχήμα \ref{fig:plot3dFunction} φαίνεται η τρισδιάστατη απεικόνιση της συνάρτησης.
-\InsertFigure{!h}{0.8}{fig:plot3dFunction}{../scripts/figures/FunctionPlot.png}{Γραφική παράσταση της f}
+\InsertFigure{!h}{0.8}{fig:plot3dFunction}{../scripts/figures/Plot_Function.png}{Γραφική παράσταση της f}
 
 Από το σχήμα μπορούμε πολύ εύκολα να διακρίνουμε ότι η συνάρτηση έχει ένα ευκρινές μέγιστο και ένα ελάχιστο στο διάστημα $x,y \in [-3, 3]$.
 Για να πάρουμε μια καλύτερη αίσθηση για το που βρίσκονται αυτά τα τοπικά ακρότατα, παρακάτω παραθέτουμε ένα γράφημα με τις ισοβαρείς καμπύλες της $f$.
-\InsertFigure{H}{0.8}{fig:plotContour}{../scripts/figures/FunctionContour.png}{Ισοβαρείς της f}
+\InsertFigure{H}{0.8}{fig:plotContour}{../scripts/figures/Plot_Contour.png}{Ισοβαρείς της f}
 
 Από το παραπάνω σχήμα \ref{fig:plotContour} φαίνεται ότι το ελάχιστο της f βρίσκεται στο αρνητικό ημιεπίπεδο των χ, κοντά στο $y = 0$
 \par
@@ -116,10 +124,10 @@
 Πριν προχωρήσουμε στα επόμενα θέματα της εργασίας και στην ανάλυση των μεθόδων υπολογισμού του ελάχιστου, θέλουμε να αναφερθούμε στις διαφορετικές τεχνικές επιλογής βήματος $\gamma_k$ και ειδικότερα για αυτή της ελαχιστοποίησης της $f(x_k + \gamma_k d_k)$ και την Armijo.
 
 \subsection{Ελαχιστοποίηση της $f(x_k + \gamma_k d_k)$}
-Η μέθοδος αυτή αναζητά την τιμή $\gamma_k$ που ελαχιστοποιεί την τιμή της συνάρτησης κατά μήκος της κατεύθυνσης $d_k$.
+Η μέθοδος αυτή αναζητά την τιμή $\gamma_k$ που ελαχιστοποιεί την τιμή της συνάρτησης κατά μήκος της κατεύθυνσης $d_k$ (αρχείο: \textbf{gamma\_minimized.m}).
 Δηλαδή, λύνουμε το μονοδιάστατο πρόβλημα:
 \boldmath
-\[\displaystyle \min{\gamma_k} f(x_k + \gamma_k d_k) \]
+\[\displaystyle \min{\gamma_k} \{f(x_k + \gamma_k d_k)\} \]
 Η κατεύθυνση $d_k$ μπορεί να είναι:
 \begin{itemize}
 	\item Η αρνητική κλίση $-\nabla f(x_k)$ (Steepest Descent).
@@ -130,127 +138,195 @@
 Στα μειονεκτήματα μπορούμε να αναφέρουμε το \textbf{υπολογιστικό κόστος}, καθώς ο υπολογισμός του $\gamma_k$ απαιτεί πολλαπλές αξιολογήσεις της $f(x)$.
 
 \subsection{Armijo rule}
-Η Armijo rule είναι μια προσαρμοστική τεχνική που επιλέγει το $\gamma_k$ για να εξασφαλίσει επαρκή μείωση της συνάρτησης.
+Η Armijo rule είναι μια προσαρμοστική τεχνική που επιλέγει το $\gamma_k$ για να εξασφαλίσει επαρκή μείωση της συνάρτησης (αρχείο: \textbf{gamma\_armijo.m}).
 Η βασική ιδέα είναι ότι η συνάρτηση πρέπει να μειώνεται "αρκετά" σε κάθε βήμα, χωρίς να χρειάζεται να υπολογίζεται το ακριβές ελάχιστο.
 Η συνθήκη του Armijo είναι:
 \boldmath
-\[ f(x_k + \gamma_k d_k) \leq f(x_k) + \sigma\gamma_k {d_k}^T*\nabla f(x_k) \]
+\[ f(x_k + \gamma_k d_k) \leq f(x_k) + \sigma\gamma_k {d_k}^T\nabla f(x_k) \]
 Όπου $\sigma \in (0, 0.1)$ είναι μια σταθερά (τυπικά $\sigma = 0.1$) και $\gamma_k$ αρχικά να ορίζεται ως 1 και να μειώνεται προοδευτικά (π.χ., $\gamma_k = \beta \cdot \gamma_k$) έως ότου ικανοποιηθεί η συνθήκη.
 \unboldmath
 \par
 Πλεονεκτήματα της μεθόδου είναι η \textbf{σταθερότητα}, καθώς αποτρέπει πολύ μεγάλα βήματα που μπορεί να αυξήσουν την τιμή της $f(x)$, αλλά και η \textbf{ανθεκτικότητα}, καθώς λειτουργεί καλά ακόμα και όταν η $f(x)$ δεν συμπεριφέρεται πολύ καλά.
 Στα μειονεκτήματα μπορούμε να αναφέρουμε την \textbf{εξάρτησή της από τις παραμέτρους} \boldmath $\sigma$ και $\beta$\unboldmath, μια κακή επιλογή των οποίων μπορεί να οδηγήσει σε αργή σύγκλιση.
 
-\section{Μέθοδος Μέγιστης Καθόδου - Steepest Descent}
-Η πρώτη μέθοδος που χρησιμοποιούμε στην εργασία (Θέμα 2), είναι η μέθοδος μέγιστης καθόδου.
+\section{Μέθοδος Μέγιστης Καθόδου - Θέμα 2}
+Η πρώτη μέθοδος που χρησιμοποιούμε στην εργασία, είναι η μέθοδος μέγιστης καθόδου (Steepest descent) (αρχείο: \textbf{method\_steepest\_descent.m}).
 Είναι μια μέθοδος πρώτης τάξης που χρησιμοποιεί την κατεύθυνση της αρνητικής κλίσης $\nabla f(x,y)$ της $f$ ως κατεύθυνση καθόδου.
 Η μέθοδος θεωρείται βασική και συχνά χρησιμοποιείται ως εισαγωγή στις μεθόδους βελτιστοποίησης.
 \par
 Η μέθοδος επιλέγει την κατεύθυνση $d_k = -\nabla f(x_k)$, η οποία είναι η κατεύθυνση της μέγιστης τοπικής μείωσης της συνάρτησης.
 Στη συνέχεια, υπολογίζεται το βήμα $\gamma_k$​ για να βρεθεί το επόμενο σημείο $x_{k+1} = x_k + \gamma_k d_k$.
 Για να χρησιμοποιήσουμε τη μέθοδο, η συνάρτηση $f$ \textbf{πρέπει να είναι συνεχής και διαφορίσιμη} και η κλίση $\nabla f$ να είναι υπολογίσιμη.
-Επίσης για την εφαρμογή της μεθόδου το αρχικό σημείο θα πρέπει να \textbf{μην είναι ακρότατο} της $f$, δηλαδή \boldmath$\nabla f(x_0) \neq 0$\unboldmath.
-\par
+Επίσης για την εφαρμογή της μεθόδου, για το αρχικό σημείο θα πρέπει \boldmath$\nabla f(x_0) \neq 0$\unboldmath.
 Όλοι οι υπολογισμοί και τα διαγράμματα για τη μέθοδο βρίσκονται στο αρχείο \textbf{Script\_2\_Steepest\_descent.m}
 
 \subsection{Σημείο εκκίνησης (0,0)}
 Για το σημείο (0, 0) η κλίση της $f$ είναι: $\nabla f(0,0) = \begin{bmatrix} 0 \\ 0 \end{bmatrix}$, με αποτέλεσμα η μέθοδος να μην μπορεί να εφαρμοστεί για κανένα τρόπο υπολογισμού βήματος.
 
 \subsection{Σημείο εκκίνησης (-1,1)}
-Για το σημείο (-1, 1) η τιμή της $f$ είναι: $f(-1, 1) = -0.135335$ και το διάνυσμα της κλίσης: $\nabla f(0,0) = \begin{bmatrix} 0.4060 \\ 0.2707 \end{bmatrix}$, επομένως μπορούμε να εφαρμόσουμε τη μέθοδο.
+Για το σημείο (-1, 1) η τιμή της $f$ είναι: $f(-1, 1) = -0.1353$ και το διάνυσμα της κλίσης: $\nabla f(0,0) = \begin{bmatrix} 0.4060 \\ 0.2707 \end{bmatrix}$, επομένως μπορούμε να εφαρμόσουμε τη μέθοδο.
+
 \par
 \underline{Σταθερό βήμα} \\
 Επιλέγοντας ακρίβεια $\epsilon = 0.0001$, εκτελούμε την μέθοδο \textit{method\_steepest\_descent()} και υπολογίζουμε τον αριθμό επαναλήψεων για διαφορετικές τιμές $\gamma_k$.
-\InsertFigure{H}{0.8}{fig:point2ItersOverGamma}{../scripts/figures/StDes_Iter_o_gamma_2.png}{Αριθμός επαναλήψεων για διαφορετικές τιμές $\gamma_k$ [Μέγιστη Κάθοδος].}
+\InsertFigure{H}{0.8}{fig:StDes_Iter_o_gamma_2}{../scripts/figures/StDes_Iter_o_gamma_2.png}{Αριθμός επαναλήψεων για διαφορετικές τιμές $\gamma_k$ [Μέγιστη Κάθοδος].}
 
-Στο παραπάνω σχήμα \ref{fig:point2ItersOverGamma} παρατηρούμε ότι για τιμές του $\gamma_k > 0.61$ η μέθοδος αποκλίνει.
+Στο παραπάνω σχήμα \ref{fig:StDes_Iter_o_gamma_2} παρατηρούμε ότι για τιμές του $\gamma_k > 0.61$ η μέθοδος αποκλίνει.
 Από την παραπάνω διαδικασία επίσης υπολογίζουμε το $\gamma_k = 0,46768$ για το οποίο η μέθοδος συγκλίνει με τα λιγότερα βήματα.
 Στο παρακάτω σχήμα \ref{fig:StDes_fixed_2} αναπαριστούμε την πορεία των σημείων καθώς συγκλίνουν στο ελάχιστο.
 \InsertFigure{H}{0.8}{fig:StDes_fixed_2}{../scripts/figures/StDes_fixed_2.png}{Σύγκλιση της μεθόδου Steepest descent [fixed $\gamma_k$].}
+
 \par
 \underline{Ελαχιστοποίηση της $f(x_k + \gamma_k d_k$)} \\
 Για την ελαχιστοποίηση της $f$, χρησιμοποιήθηκε η bisection από την προηγούμενη εργασία, η οποία τροποποιήθηκε ώστε δέχεται functions και όχι symbolic expressions.
 \InsertFigure{H}{0.8}{fig:StDes_minimized_2}{../scripts/figures/StDes_minimized_2.png}{Σύγκλιση της μεθόδου Steepest descent [minimized f].}
 Από το γράφημα φαίνεται τόσο ότι η μέθοδος συγκλίνει κοντά στο ελάχιστο γρηγορότερα, όσο και ότι πραγματοποιεί “διορθώσεις πορείας”.
+
 \par
 \underline{Armijo rule} \\
 Για τη μέθοδο η βασική ιδέα είναι να ξεκινήσει ο αλγόριθμος από ένα μεγάλο $\gamma_k = 1$ και συνεχώς να μειώνεται με βάση τον κανόνα Armijo.
 Μετά από ένα tuning των παραμέτρων της μεθόδου καταλήξαμε στα $\beta=0.4, \sigma=0.1$
 \InsertFigure{H}{0.8}{fig:StDes_armijo_2}{../scripts/figures/StDes_armijo_2.png}{Σύγκλιση της μεθόδου Steepest descent [armijo rule].}
 
+\par
+\underline{Σύγκριση} \\
+Στο παρακάτω διάγραμμα φαίνονται τα βήματα αλλά ο ρυθμός σύγκλισης της κάθε μεθόδου.
+\InsertFigure{H}{0.8}{fig:StDes_compare_2}{../scripts/figures/StDes_compare_2.png}{Ευκλείδεια απόσταση από το ελάχιστο για κάθε μέθοδο υπολογισμού $\gamma_k$ [Steepest descent].}
+
 \subsection{Σημείο εκκίνησης (1,-1)}
-Για το σημείο (1, -1) η τιμή της $f$ είναι: $f(1, -1) = -0.135335$ και το διάνυσμα της κλίσης: $\nabla f(0,0) = \begin{bmatrix} 0.4060 \\ 0.2707 \end{bmatrix}$, επομένως μπορούμε να εφαρμόσουμε τη μέθοδο.
+Για το σημείο (1, -1) η τιμή της $f$ είναι: $f(1, -1) = -0.1353$ και το διάνυσμα της κλίσης: $\nabla f(0,0) = \begin{bmatrix} 0.4060 \\ 0.2707 \end{bmatrix}$, επομένως μπορούμε να εφαρμόσουμε τη μέθοδο.
+
 \par
 \underline{Σταθερό βήμα} \\
-Για σταθερό βήμα εκτελέσαμε διαδοχικά τη μέθοδο \textit{method\_steepest\_descent()} για να υπολογίσουμε τον αριθμό επαναλήψεων για διαφορετικές τιμές $\gamma_k$, όμως σε καμία τιμή ο αλγόριθμος δεν συγκλίνει.
-Ακόμα και για μεγάλες τιμές του $\gamma_k$, ο αλγόριθμος εγκλωβίζεται στο δεξιό ημιεπίπεδο.
+Για σταθερό βήμα εκτελέσαμε διαδοχικά τη μέθοδο \textit{method\_steepest\_descent()} για να υπολογίσουμε τον αριθμό επαναλήψεων για διαφορετικές τιμές $\gamma_k$, όμως σε καμία τιμή ο αλγόριθμος δεν καταφέρνει να συγκλίνει.
+Ακόμα και για μεγάλες τιμές βήματος, όπως για παράδειγμα $\gamma_k = 1$, ο αλγόριθμος εγκλωβίζεται στο δεξιό ημιεπίπεδο, όπως φαίνεται και στο διάγραμμα \ref{fig:StDes_fixed_3} παρακάτω.
 \InsertFigure{H}{0.8}{fig:StDes_fixed_3}{../scripts/figures/StDes_fixed_3.png}{Μη σύγκλιση της μεθόδου Steepest descent [Fixed step].}
+
 \par
 \underline{Ελαχιστοποίηση της $f(x_k + \gamma_k d_k$)} \\
 \InsertFigure{H}{0.8}{fig:StDes_minimized_3}{../scripts/figures/StDes_minimized_3.png}{Σύγκλιση της μεθόδου Steepest descent [minimized f].}
-Από το γράφημα φαίνεται ότι η μέθοδος συγκλίνει, καταφέρνοντας να περάσει την περιοχή με μηδενικές κλίσεις κοντά στον άξονα των $\psi$.
+Για υπολογισμό βήματος ελαχιστοποιώντας τη $f(x_k + \gamma_k d_k)$, όπως φαίνεται και από το γράφημα \ref{fig:StDes_minimized_3}, η μέθοδος συγκλίνει, καταφέρνοντας να περάσει την περιοχή με μηδενικές κλίσεις κοντά στον άξονα των y.
+
 \par
 \underline{Armijo rule} \\
 \InsertFigure{H}{0.8}{fig:StDes_armijo_3}{../scripts/figures/StDes_armijo_3.png}{Μη σύγκλιση της μεθόδου Steepest descent [armijo rule].}
-Αντίθετα η μέθοδος armijo δεν συγκλίνει, καθώς εγκλωβίζεται στο δεξιό ημιεπίπεδο.
+Αντίθετα η μέθοδος armijo δεν συγκλίνει, καθώς και αυτή εγκλωβίζεται στο δεξιό ημιεπίπεδο.
 
-\section{Μέθοδος Newton}
-Η δεύτερη μέθοδος που χρησιμοποιούμε στην εργασία (Θέμα 3), είναι η μέθοδος Newton.
-Η μέθοδος χρησιμοποιεί πληροφορίες δεύτερης τάξης (Hessian) για τη βελτίωση της κατεύθυνσης καθόδου.
-Αν η συνάρτηση είναι τετραγωνική, η μέθοδος συγκλίνει σε ένα βήμα.
+\section{Μέθοδος Newton - Θέμα 3}
+Η δεύτερη μέθοδος που χρησιμοποιούμε στην εργασία, είναι η μέθοδος Newton (αρχείο: \textbf{method\_newton.m}).
+Η μέθοδος χρησιμοποιεί πληροφορίες δεύτερης τάξης (Εσσιανό - Hessian) για τη βελτίωση της κατεύθυνσης καθόδου.
 Η μέθοδος ορίζει την κατεύθυνση
 \boldmath\[d_k = -{H_k}^{-1}\nabla f(x_k)\]\unboldmath
 Όπου $H_k$ είναι ο Εσσιανός πίνακας της $f$ στο $x_k$.
 Το επόμενο σημείο υπολογίζεται ως
 \boldmath\[x_{k+1} = x_k + \gamma_k d_k\]\unboldmath
-Με κατάλληλο υπολογισμό του βήματος.
-Για να λειτουργήσει η μέθοδος η $f$ πρέπει να είναι \textbf{δύο φορές διαφορίσιμη} και ο Εσσιανός \boldmath$H_k$\unboldmath να είναι \textbf{θετικά ορισμένος και αντιστρέψιμος}.
+Ο υπολογισμός βήματος γίνεται παλι είτε με σταθερό βήμα, είτε με ελαχιστοποίηση της $f(x_k + \gamma_k d_k)$, είτε με τον κανόνα armijo. 
+Για να λειτουργήσει η μέθοδος, η $f$ πρέπει να είναι \textbf{δύο φορές διαφορίσιμη} και ο Εσσιανός \boldmath$H_k$\unboldmath να είναι \textbf{θετικά ορισμένος και αντιστρέψιμος}.
 \par
 Στα πλεονεκτήματα της μεθόδου είναι η \textbf{ταχύτερη σύγκλιση} από την Steepest Descent για κυρτές συναρτήσεις και το γεγονός ότι εκμεταλλεύεται την \textbf{πληροφορία καμπυλότητας} της συνάρτησης.
 Όμως είναι υπολογιστικά δαπανηρή και δεν είναι ανθεκτική σε μη κυρτές συναρτήσεις ή σε περιπτώσεις όπου ο Εσσιανός είναι κακώς ορισμένος.
-
 Όλοι οι υπολογισμοί για τη μέθοδο βρίσκονται στο αρχείο \textbf{Script\_3\_Newton.m}
 
 \subsection{Σημείο εκκίνησης (0,0)}
 Για το σημείο $x_k = (0, 0)$ η κλίση της $f$ είναι: $\nabla f(x_k) = \begin{bmatrix} 0 \\ 0 \end{bmatrix}$ και ο εσσιανός $H(x_k) = \begin{bmatrix} 0 & 0 \\ 0 & 0 \end{bmatrix}$ με αποτέλεσμα η μέθοδος και εδώ να μην μπορεί να εφαρμοστεί για κανένα τρόπο υπολογισμού βήματος.
 
 \subsection{Σημείο εκκίνησης (-1,1)}
-Για το σημείο $x_k = (-1, 1)$ η κλίση της $f$ είναι: $\nabla f(x_k) = \begin{bmatrix} 0.406 \\ 0.270 \end{bmatrix}$ και ο εσσιανός $H(x_k) = \begin{bmatrix} -0.270 & -0.812 \\ -0.812 & -0.270 \end{bmatrix}$ με ιδιοτιμές $\lambda = \begin{bmatrix} -1.082 \\ 0.541 \end{bmatrix}$.
+Για το σημείο $x_k = (-1, 1)$ η κλίση της $f$ είναι: $\nabla f(x_k) = \begin{bmatrix} 0.4060 \\ 0.2707 \end{bmatrix}$ και ο εσσιανός $H(x_k) = \begin{bmatrix} -0.2707 & -0.8120 \\ -0.8120 & -0.2707 \end{bmatrix}$ με ιδιοτιμές $\lambda = \begin{bmatrix} -1.0827 \\ 0.5413 \end{bmatrix}$.
 Από τα παραπάνω προκύπτει πως ο Εσσιανός είναι αόριστος και άρα δεν μπορεί να εφαρμοστεί η μέθοδος, για κανένα τρόπο υπολογισμού βήματος.
 
 \subsection{Σημείο εκκίνησης (1,-1)}
-Για το σημείο $x_k = (1, -1)$ η κλίση της $f$ είναι: $\nabla f(x_k) = \begin{bmatrix} 0.406 \\ 0.270 \end{bmatrix}$ και ο εσσιανός $H(x_k) = \begin{bmatrix} 0.270 & 0.812 \\ 0.812 & 0.270 \end{bmatrix}$ με ιδιοτιμές $\lambda = \begin{bmatrix} -0.541 \\ 1.082 \end{bmatrix}$.
+Για το σημείο $x_k = (1, -1)$ η κλίση της $f$ είναι: $\nabla f(x_k) = \begin{bmatrix} 0.4060 \\ 0.2707 \end{bmatrix}$ και ο εσσιανός $H(x_k) = \begin{bmatrix} 0.2707 & 0.8120 \\ 0.8120 & 0.2707 \end{bmatrix}$ με ιδιοτιμές $\lambda = \begin{bmatrix} -0.5413 \\ 1.0827 \end{bmatrix}$.
 Και εδώ, από τα παραπάνω προκύπτει πως ο Εσσιανός είναι αόριστος και άρα δεν μπορεί να εφαρμοστεί η μέθοδος, για κανένα τρόπο υπολογισμού βήματος.
 
-\section{Μέθοδος Levenberg-Marquardt}
-Η τελευταία μέθοδος που χρησιμοποιούμε στην εργασία (Θέμα 4), είναι η μέθοδος Levenberg-Marquardt.
+\section{Μέθοδος Levenberg-Marquardt - Θέμα 4}
+Η τελευταία μέθοδος που χρησιμοποιούμε στην εργασία, είναι η μέθοδος Levenberg-Marquardt (αρχείο: \textbf{method\_lev\_mar.m}).
 Πρόκειται για μια τροποποιημένη έκδοση της μεθόδου Newton, η οποία εισάγει έναν παράγοντα απόσβεσης για τη σταθεροποίηση όταν ο εσσιανός δεν είναι θετικά ορισμένος.
-Για το λόγο αυτό χρησιμοποιεί ένας προσαρμοσμένος εσσιανός $H_k' = H_k + \mu_k I$, όπου $\mu_k > 0$ ένας παράγοντας απόσβεσης.
-Για μεγάλες τιμές του $\mu_k$ η μέθοδος συμπεριφέρεται σαν Gradient Descent.
-
-%Απαιτήσεις:
-%Η f(x)f(x) πρέπει να είναι δύο φορές διαφορίσιμη.
-%Υπολογισμός του λλ απαιτεί προσεκτική επιλογή παραμέτρων.
-%Περιορισμοί:
-%Η απόδοση εξαρτάται από την επιλογή του αρχικού λλ.
-%Μπορεί να παρουσιάσει αργή σύγκλιση σε προβλήματα χωρίς κυρτότητα.
-%Πλεονεκτήματα:
-%Σταθερή ακόμη και για κακώς ορισμένα Hessians.
-%Λειτουργεί καλά σε προβλήματα που συνδυάζουν γραμμικές και μη γραμμικές εξαρτήσεις.
-%Μειονεκτήματα:
-%Υψηλότερο υπολογιστικό κόστος σε σχέση με το Steepest Descent.
-
+\boldmath
+Για το λόγο αυτό χρησιμοποιείται ένας προσαρμοσμένος εσσιανός $H_k' = H_k + \mu_k I$, όπου $\mu_k > 0$ ένας παράγοντας, τέτοιος ώστε ο $H_k'$ να είναι θετικά ορισμένος.
+Έτσι σε κάθε βήμα, γίνεται έλεγχος των ιδιοτιμών του εσσιανού.
+Αν υπάρχει ιδιοτιμή μικρότερη ή ίση με το 0, τότε υπολογίζεται ο παράγοντας \[\mu_k = \abs{\min\{\lambda(H_k)\}} + \epsilon \]
+Όπου $\epsilon > 0$ ένας παράγοντας προσαύξησης ώστε να βεβαιωθούμε πως ο παραγόμενος προσαρμοσμένος Εσσιανός $H_k' = H_k + \mu_k I$ είναι θετικά ορισμένος και δεν έχει κάποια μηδενική ιδιοτιμή.
+\unboldmath
+Να σημειώσουμε εδώ πως για μεγάλες τιμές του $\mu_k$ η μέθοδος συμπεριφέρεται σαν Gradient Descent, ενώ για μικρές προσεγγίζει τη μέθοδο Newton.
+\par
+Για να λειτουργήσει η μέθοδος, η $f$ πρέπει να είναι \textbf{δύο φορές διαφορίσιμη}.
+Επίσης σημαντικό ρόλο στην ευστάθεια της μεθόδου παίζει η επιλογή το $\epsilon$ και άρα κατ' επέκταση και του $\mu_k$.
+Στα πλεονεκτήματα της μεθόδου έχουμε την σταθερότητα ακόμα και για κακώς ορισμένους Εσσιανούς πίνακες, η οποία όμως έρχεται με υπολογιστικό κόστος.
 Όλοι οι υπολογισμοί για τη μέθοδο βρίσκονται στο αρχείο \textbf{Script\_4\_LevMar.m}
 
-\section{Σύγκριση των μεθόδων}
-Εκτελώντας όλους του αλγόριθμους για τα ίδια δεδομένα, \textbf{για τον αριθμό επαναλήψεων} έχουμε: \\
+\subsection{Σημείο εκκίνησης (0,0)}
+Για το σημείο $x_k = (0, 0)$ η κλίση της $f$ είναι: $\nabla f(x_k) = \begin{bmatrix} 0 \\ 0 \end{bmatrix}$ και ο εσσιανός $H(x_k) = \begin{bmatrix} 0 & 0 \\ 0 & 0 \end{bmatrix}$ με αποτέλεσμα η μέθοδος και εδώ να μην μπορεί να εφαρμοστεί για κανένα τρόπο υπολογισμού βήματος.
 
-\begin{itemize}
-    \item ...
-\end{itemize}
+\subsection{Σημείο εκκίνησης (-1,1)}
+Για το σημείο $x_k = (0, 0)$ η κλίση της $f$ είναι: $\nabla f(x_k) = \begin{bmatrix} 0.4060 \\ 0.2707 \end{bmatrix}$ και η μέθοδος μπορεί να εφαρμοστεί.
+\par
+\underline{Σταθερό βήμα} \\
+Επιλέγοντας ακρίβεια $\epsilon = 0.0001$, προσπαθούμε να βρούμε τον συντελεστή $\epsilon$ για τον οποίο η μέθοδος είναι σταθερή.
+Κατόπιν εκτελούμε την μέθοδο \textit{method\_lev\_mar()} και υπολογίζουμε τον αριθμό επαναλήψεων για διαφορετικές τιμές $\gamma_k$.
+Μετά από tuning, καταλήξαμε στην τιμή $\epsilon = 0.3$ για την οποία αναζητήσαμε και υπολογίσαμε το βήμα $\gamma_k = 1.4152$ με τη γρηγορότερη σύγκλιση.
+\InsertFigure{H}{0.8}{fig:LevMar_Iter_o_gamma_2}{../scripts/figures/LevMar_Iter_o_gamma_2.png}{Αριθμός επαναλήψεων για διαφορετικές τιμές $\gamma_k$ [Μέθοδος L-M].}
 
+Στο παρακάτω σχήμα \ref{fig:LevMar_fixed_2} αναπαριστούμε την πορεία των σημείων καθώς συγκλίνουν στο ελάχιστο.
+Παρατηρούμε ωστόσο ότι η μέθοδος κάνει αρκετή "διόρθωση πορείας" μέχρι τη σύγκλιση.
+\InsertFigure{H}{0.8}{fig:LevMar_fixed_2}{../scripts/figures/LevMar_fixed_2.png}{Σύγκλιση της μεθόδου Levenberg-Marquardt [fixed $\gamma_k$].}
 
-\section{Συμπεράσματα}
-Οι μέθοδοι της παρούσας εργασίας αποτελούν βασικές τεχνικές για την εύρεση του τοπικού ελαχίστου ...
+\par
+\underline{Ελαχιστοποίηση της $f(x_k + \gamma_k d_k$)} \\
+Για την ελαχιστοποίηση της $f$, χρησιμοποιήθηκε και εδώ η bisection από την προηγούμενη εργασία.
+\InsertFigure{H}{0.8}{fig:LevMar_minimized_2}{../scripts/figures/LevMar_minimized_2.png}{Σύγκλιση της μεθόδου Levenberg-Marquardt [minimized f].}
+Από το γράφημα φαίνεται τόσο ότι η μέθοδος συγκλίνει κοντά στο ελάχιστο, όσο και ότι πραγματοποιεί “διορθώσεις πορείας”, μικρότερες όμως από αυτές του σταθερού βήματος.
+
+\par
+\underline{Armijo rule} \\
+Για τη μέθοδο η βασική ιδέα και εδώ είναι να ξεκινήσει ο αλγόριθμος από ένα μεγάλο $\gamma_k = 1$ και συνεχώς να μειώνεται με βάση τον κανόνα Armijo.
+Μετά από ένα tuning των παραμέτρων της μεθόδου καταλήξαμε στα $\beta=0.4, \sigma=0.1$
+\InsertFigure{H}{0.8}{fig:LevMar_armijo_2}{../scripts/figures/LevMar_armijo_2.png}{Σύγκλιση της μεθόδου Levenberg-Marquardt [armijo rule].}
+
+\par
+\underline{Σύγκριση} \\
+Στο παρακάτω διάγραμμα φαίνονται τα βήματα αλλά ο ρυθμός σύγκλισης της κάθε μεθόδου.
+\InsertFigure{H}{0.8}{fig:LevMar_compare_2}{../scripts/figures/LevMar_compare_2.png}{Ευκλείδεια απόσταση από το ελάχιστο για κάθε μέθοδο υπολογισμού $\gamma_k$ [Levenberg-Marquardt].}
+
+\subsection{Σημείο εκκίνησης (1,-1)}
+Για το σημείο (1, -1) η τιμή της $f$ είναι: $f(1, -1) = -0.1353$ και το διάνυσμα της κλίσης: $\nabla f(0,0) = \begin{bmatrix} 0.4060 \\ 0.2707 \end{bmatrix}$, επομένως μπορούμε να εφαρμόσουμε τη μέθοδο.
+
+\par
+\underline{Σταθερό βήμα} \\
+Για σταθερό βήμα εκτελέσαμε διαδοχικά τη μέθοδο \textit{method\_lev\_mar()} για να υπολογίσουμε τον αριθμό επαναλήψεων για διαφορετικές τιμές $\gamma_k$, όμως σε καμία τιμή ο αλγόριθμος δεν καταφέρνει να συγκλίνει.
+Παρακάτω (σχήμα \ref{fig:LevMar_fixed_3}), παραθέτουμε ένα παράδειγμα $\gamma_k = 0.1$, όπου ο αλγόριθμος εγκλωβίζεται στο δεξιό ημιεπίπεδο.
+\InsertFigure{H}{0.8}{fig:LevMar_fixed_3}{../scripts/figures/LevMar_fixed_3.png}{Μη σύγκλιση της μεθόδου Levenberg-Marquardt [Fixed step].}
+
+\par
+\underline{Ελαχιστοποίηση της $f(x_k + \gamma_k d_k$)} \\
+\InsertFigure{H}{0.8}{fig:LevMar_minimized_3}{../scripts/figures/LevMar_minimized_3.png}{Σύγκλιση της μεθόδου Levenberg-Marquardt [minimized f].}
+Για υπολογισμό βήματος ελαχιστοποιώντας τη $f(x_k + \gamma_k d_k)$, όπως φαίνεται και από το γράφημα \ref{fig:LevMar_minimized_3}, η μέθοδος συγκλίνει, καταφέρνοντας να περάσει την περιοχή με μηδενικές κλίσεις κοντά στον άξονα των $\psi$.
+
+\par
+\underline{Armijo rule} \\
+\InsertFigure{H}{0.8}{fig:LevMar_armijo_3}{../scripts/figures/LevMar_armijo_3.png}{Μη σύγκλιση της μεθόδου Levenberg-Marquardt [armijo rule].}
+Αντίθετα η μέθοδος armijo δεν συγκλίνει, καθώς και αυτή εγκλωβίζεται στο δεξιό ημιεπίπεδο.
+
+\section{Σύγκριση των μεθόδων - Συμπεράσματα}
+Εκτελώντας όλους του αλγόριθμους σε συνδυασμό με όλες τις μεθόδους, για το σημείο (-1, 1), \textbf{για τον αριθμό επαναλήψεων} έχουμε: \\
+
+\noindent
+\renewcommand{\arraystretch}{1.2}
+\begin{tabularx}{\textwidth}{%
+		>{\raggedleft\arraybackslash}m{0.21\textwidth} |
+		>{\centering\arraybackslash}m{0.21\textwidth} |
+		>{\centering\arraybackslash}m{0.21\textwidth} |
+		>{\centering\arraybackslash}m{0.21\textwidth}
+	}
+	Μέθοδος & Σταθερό βήμα & Ελαχιστοποίηση της f & Armijo \\
+	\hline
+	St. Descent & 11 & 10 & 11 \\
+	Newton      &  - &  - &  - \\
+	Lev-Mar     &  7 &  7 &  9
+\end{tabularx} \\ [2ex]
+Από τον παραπάνω πίνακα, αλλά και από τη πορεία σύγκλισης συμπεραίνουμε πως ο αλγόριθμος Lev-Mar είναι ο πιο αποδοτικός.
+Επίσης σε συνδυασμό με τον τρόπο επιλογής βήματος το οποίο ελαχιστοποιεί τη $f(x_k + \gamma_kd_k)$ ο αλγόριθμος εκτός από γρηγορότερος, καταφέρνει να απεγκλωβίζεται από περιοχές με πολύ μικρή κλίση.
+Αντίθετα οι μέθοδοι με σταθερό βήμα και armijo, φαίνεται ότι παρουσιάζουν τοπικό χαρακτήρα, καθώς εγκλωβίζονται εύκολα, με αποτέλεσμα να μην υπάρχει εγγύηση ότι το σημείο που βρέθηκε είναι ολικό ελάχιστο.
+Ακόμα η μέθοδος Newton δεν φάνηκε να μπορεί να χρησιμοποιηθεί στην παρούσα εργασία, καθώς η επιλογή των σημείων ήταν τέτοια που είτε η κλίση ήταν 0, είτε ο Εσσιανός κακώς ορισμένος.
+Αυτό φυσικά αναδεικνύει απλώς την μικρότερη εφαρμοσιμότητα του αλγορίθμου.
 
 \end{document}

Work2/scripts/GivenEnv.m

@@ -15,6 +15,10 @@ fun         = matlabFunction(fexpr, 'Vars', [x, y]);        % Function
 grad_fun    = matlabFunction(grad_fexpr, 'Vars', [x, y]); % Gradient
 hessian_fun = matlabFunction(hessian_fexpr, 'Vars', [x, y]);  % Hessian
 
+% Minimum reference
+Freference = @(x) x(1).^5 .* exp(-x(1).^2 - x(2).^2);
+[Xmin, Fmin] = fminsearch(Freference, [-1, -1]);
+
 % Amijo globals
 global amijo_beta; % Step reduction factor in [0.1, 0.5] (typical range: [0.1, 0.8])
 global amijo_sigma; % Sufficient decrease constant in [1e-5, 0.1] (typical range: [0.01, 0.3])

Work2/scripts/Script_1_Plots.m

@@ -8,7 +8,7 @@ GivenEnv
 %
 
 % 3d plot the function
-plot3dFun(fun, [-3, 3], [-3, 3], 100, title_fun, "figures/FunctionPlot.png");
+plot3dFun(fun, [-3, 3], [-3, 3], 100, title_fun, "figures/Plot_Function.png");
 
 % Plot isobaric lines
-plotContour(fun, [-3, 3], [-3, 3], 100, title_fun, "figures/FunctionContour.png");
+plotContour(fun, [-3, 3], [-3, 3], 100, title_fun, "figures/Plot_Contour.png");

Work2/scripts/Script_2_Steepest_descent.m

@@ -63,6 +63,9 @@ fprintf('Armijo step:    Initial point (%f, %f), steps:%d, Final (x,y)=(%f, %f),
 plotPointsOverContour(x_armijo, fun, [-2, 0], [-2, 2], 100, point_str + ": Steepest descent Armijo method", "figures/StDes_armijo_" + point + ".png");  
 disp(' ');
 
+% Compare methods
+plotConvCompare(x_fixed, "Fixed", x_minimized, "Minimized", x_armijo, "Armijo", Xmin, "Convergence compare", "figures/StDes_compare_" + point + ".png");
+
 % Point x0 = (1, -1)
 % =========================================================================
 point = 3;
@@ -111,4 +114,3 @@ fprintf('Armijo step:    Initial point (%f, %f), steps:%d, Final (x,y)=(%f, %f),
 plotPointsOverContour(x_armijo, fun, [-1, 2], [-2, 2], 100, point_str + ": Steepest descent Armijo method", "figures/StDes_armijo_" + point + ".png");
 
 
-

Work2/scripts/Script_3_Newton.m

@@ -40,4 +40,3 @@ ev = eig(hf);
 fprintf('Initial point (%d, %d), f = %f, grad = [%f;%f], hessian = [%f %f ; %f %f]. Eigenvalues= [%f, %f], Can NOT use method\n', x0, f, gf, hf, ev);
 
 
-

Work2/scripts/Script_4_LevMar.m

@@ -43,6 +43,7 @@ for g = n
     end
     j = j + 1;
 end
+plotItersOverGamma(n, k, "Iteration for different $\gamma$ values", "figures/LevMar_Iter_o_gamma_" + point + ".png");
 
 [~, j] = min(k);
 gamma_fixed_step = n(j);
@@ -51,9 +52,9 @@ gamma_fixed_step = n(j);
 fprintf('Fixed step:     Initial point (%f, %f), steps:%d, Final (x,y)=(%f, %f), f(x,y)=%f\n', x0, kk, x_fixed(end, :), f_fixed(end));
 plotPointsOverContour(x_fixed, fun, [-3, 0], [-2, 2], 100, point_str + ": Levenberg-Marquardt $\gamma$ = " + gamma_fixed_step, "figures/LevMar_fixed_" + point + ".png");
 
-[x_fixed, f_fixed, kk] = method_lev_mar(fun, grad_fun, hessian_fun, 0.3, x0, tol, max_iter, 'minimized');
-fprintf('Minimized f(g): Initial point (%f, %f), steps:%d, Final (x,y)=(%f, %f), f(x,y)=%f\n', x0, kk, x_fixed(end, :), f_fixed(end));
-plotPointsOverContour(x_fixed, fun, [-3, 0], [-2, 2], 100, point_str + ": Levenberg-Marquardt minimized $f(x_k + \gamma_kd_k)$", "figures/LevMar_minimized_" + point + ".png");
+[x_minimized, f_minimized, kk] = method_lev_mar(fun, grad_fun, hessian_fun, 0.3, x0, tol, max_iter, 'minimized');
+fprintf('Minimized f(g): Initial point (%f, %f), steps:%d, Final (x,y)=(%f, %f), f(x,y)=%f\n', x0, kk, x_minimized(end, :), f_minimized(end));
+plotPointsOverContour(x_minimized, fun, [-3, 0], [-2, 2], 100, point_str + ": Levenberg-Marquardt minimized $f(x_k + \gamma_kd_k)$", "figures/LevMar_minimized_" + point + ".png");
 
 % Armijo Rule
 
@@ -63,9 +64,11 @@ amijo_sigma = 0.1; % typical range: [0.01, 0.3]
 
 [x_armijo, f_armijo, kk] = method_lev_mar(fun, grad_fun, hessian_fun, 0.3, x0, tol, max_iter, 'armijo');
 fprintf('Armijo step:    Initial point (%f, %f), steps:%d, Final (x,y)=(%f, %f), f(x,y)=%f\n', x0, kk, x_armijo(end, :), f_armijo(end));
-plotPointsOverContour(x_armijo, fun, [-3, 0], [-2, 2], 100, point_str + ": Levenberg-Marquardt Armijo method", "figures/StDes_armijo_" + point + ".png");  
+plotPointsOverContour(x_armijo, fun, [-3, 0], [-2, 2], 100, point_str + ": Levenberg-Marquardt Armijo method", "figures/LevMar_armijo_" + point + ".png");  
 disp(' ');
 
+% Compare methods
+plotConvCompare(x_fixed, "Fixed", x_minimized, "Minimized", x_armijo, "Armijo", Xmin, "Convergence compare", "figures/LevMar_compare_" + point + ".png");
 
 % Point x0 = (1, -1)
 % =========================================================================
@@ -112,5 +115,5 @@ amijo_sigma = 0.1; % typical range: [0.01, 0.3]
 
 [x_armijo, f_armijo, kk] = method_lev_mar(fun, grad_fun, hessian_fun, 0.3, x0, tol, max_iter, 'armijo');
 fprintf('Armijo step:    Initial point (%f, %f), steps:%d, Final (x,y)=(%f, %f), f(x,y)=%f\n', x0, kk, x_armijo(end, :), f_armijo(end));
-plotPointsOverContour(x_armijo, fun, [-3, 2], [-2, 2], 100, point_str + ": Levenberg-Marquardt Armijo method", "figures/StDes_armijo_" + point + ".png");  
-disp(' ');
+plotPointsOverContour(x_armijo, fun, [-3, 2], [-2, 2], 100, point_str + ": Levenberg-Marquardt Armijo method", "figures/LevMar_armijo_" + point + ".png");  
+disp(' ');

Work2/scripts/figures/ReadME.md

@@ -1 +0,0 @@
-File to keep directory

Work2/scripts/fmin_bisection.m

@@ -46,4 +46,4 @@ while b(k) - a(k) > lambda
     end
 end
 
-end
+end

Work2/scripts/gamma_armijo.m

@@ -30,4 +30,4 @@ function [gamma] = gamma_armijo(f, grad_f, dk, xk)
         end
     end
 
-end
+end

Work2/scripts/gamma_minimized.m

@@ -21,4 +21,4 @@ function [gamma] = gamma_minimized(f, ~, dk, xk)
 
     % find g that minimizes fmin
     %gamma = fminbnd(g, 0, 1);
-end
+end

Work2/scripts/method_lev_mar.m

@@ -2,7 +2,9 @@ function [x_vals, f_vals, k] = method_lev_mar(f, grad_f, hessian_f, e, xk, tol,
     % f:         Objective function
     % grad_f:    Gradient of the function
     % hessian_f: Hessian of the function
-    % e:         mu offset for hessian damping H' = H_k + mI
+    % e:         Offset for hessian damping Hk' = Hk + mI
+    %            - when:  Hk not positive defined
+    %            - Where: m = abs(min(eig(Hk))) + e
     % xk:        Initial point [xk, yk]
     % tol:       Tolerance for stopping criterion
     % max_iter:  Maximum number of iterations
@@ -35,10 +37,10 @@ function [x_vals, f_vals, k] = method_lev_mar(f, grad_f, hessian_f, e, xk, tol,
         % Check if hessian is not positive defined
         lmin = min(eig(hess));
         if lmin <= 0
+            % Select m with offset to stear hess to positive eigenvalues
             m = abs(lmin) + e;
             mI = m * eye(size(hess));
-            nev = eig(hess + mI);
-            if min(nev) <= 0
+            if min(eig(hess + mI)) <= 0 % Fail-check
                 warning('Can not normalize hessian matrix.');
             end
         end
@@ -56,4 +58,4 @@ function [x_vals, f_vals, k] = method_lev_mar(f, grad_f, hessian_f, e, xk, tol,
         x_vals = [x_vals; x_next]; % Store values
         f_vals = [f_vals; f_next]; % Store function values
     end
-end
+end

Work2/scripts/method_newton.m

@@ -44,4 +44,4 @@ function [x_vals, f_vals, k] = method_newton(f, grad_f, hessian_f, xk, tol, max_
         x_vals = [x_vals; x_next]; % Store values
         f_vals = [f_vals; f_next]; % Store function values
     end
-end
+end

Work2/scripts/method_steepest_descent.m

@@ -40,4 +40,4 @@ function [x_vals, f_vals, k] = method_steepest_descent(f, grad_f, xk, tol, max_i
         x_vals = [x_vals; x_next]; % Store values
         f_vals = [f_vals; f_next]; % Store function values
     end
-end
+end

Work2/scripts/plotConvCompare.m

@@ -0,0 +1,54 @@
+function plotConvCompare(points_1, title_1, points_2, title_2, points_3, title_3, Min_point, plot_title, filename)
+    % 3D plots a function
+    % points:       The points to plot
+    % contur_fun:   The function for contour plot
+    % x_lim:        The range for x axis. ex: [-2, 2]
+    % y_lim:        The range for y axis. ex: [0, 2]
+    % size:         The number of points for each axis
+    % plot_title:   The latex title for the plot
+    % filename:     The filename to save the plot (if exists)
+    %
+    
+    global image_width, 
+    global image_height;
+
+    distances_1 = sqrt((points_1(:,1) - Min_point(1)).^2 + (points_1(:,2) - Min_point(2)).^2);
+    distances_2 = sqrt((points_2(:,1) - Min_point(1)).^2 + (points_2(:,2) - Min_point(2)).^2);
+    distances_3 = sqrt((points_3(:,1) - Min_point(1)).^2 + (points_3(:,2) - Min_point(2)).^2);
+
+    
+    % 2D plot
+    figure('Name', 'Convergence compare', 'NumberTitle', 'off');
+    set(gcf, 'Position', [100, 100, image_width, image_height]); % Set the figure size
+    title(plot_title, 'Interpreter', 'latex', 'FontSize', 16); % Title of the plot
+    
+    % One
+    subplot(3, 1, 1);
+    plot(distances_1, '-o');
+    % Customize the plot
+    ylabel(title_1, 'Interpreter', 'none');
+    xlabel('Step');
+    grid on
+
+    % One
+    subplot(3, 1, 2);
+    plot(distances_2, '-o');
+    % Customize the plot
+    ylabel(title_2, 'Interpreter', 'none');
+    xlabel('Step');
+    grid on
+
+    % One
+    subplot(3, 1, 3);
+    plot(distances_3, '-o');
+    % Customize the plot
+    ylabel(title_3, 'Interpreter', 'none');
+    xlabel('Step');
+    grid on
+
+
+    % save the figure
+    if strcmp(filename, '') == 0
+        print(gcf, filename, '-dpng', '-r300');
+    end
+end