Compare commits
10
Commits
v1.0
...
be42c39adc
| Author | SHA1 | Date | |
|---|---|---|---|
 hoo2
|
be42c39adc | ||
|
hoo2
|
65363b3717 | ||
|
hoo2
|
ed9807b1ca | ||
|
hoo2
|
d5c8ed5e68 | ||
 Christos Choutouridis
|
9565ac2538 | ||
|
Christos Choutouridis
|
2b2df13940 | ||
|
hoo2
|
3b40c20ec0 | ||
|
Christos Choutouridis
|
65fd57a70f | ||
|
hoo2
|
c50bb6e03a | ||
|
hoo2
|
9f8702729d |
+1
-1
@@ -6,4 +6,4 @@
|
||||
|
||||
# Matlab related
|
||||
*.m~
|
||||
|
||||
*.asv
|
||||
|
||||
Binary file not shown.
@@ -0,0 +1,311 @@
|
||||
%
|
||||
% Optimization Techniques Work 3 report
|
||||
%
|
||||
% authors:
|
||||
% Χρήστος Χουτουρίδης ΑΕΜ 8997
|
||||
% cchoutou@ece.auth.gr
|
||||
|
||||
|
||||
\documentclass[a4paper, 11pt]{AUTHReport}
|
||||
|
||||
% Document configuration
|
||||
\AuthorName{Χρήστος Χουτουρίδης}
|
||||
\AuthorAEM{8997}
|
||||
\AuthorMail{cchoutou@ece.auth.gr}
|
||||
|
||||
%\CoAuthorName{CoAuthor Name}
|
||||
%\CoAuthorAEM{AEM}
|
||||
%\CoAuthorMail{CoAuthor Mail}
|
||||
|
||||
% \WorkGroup{Ομάδα Χ}
|
||||
|
||||
\DocTitle{3η Εργαστηριακή Άσκηση}
|
||||
\DocSubTitle{Μέθοδος Μέγιστης Καθόδου με Προβολή}
|
||||
|
||||
\Department{Τμήμα ΗΜΜΥ. Τομέας Ηλεκτρονικής}
|
||||
\ClassName{Τεχνικές Βελτιστοποίησης}
|
||||
|
||||
\InstructorName{Γ. Ροβιθάκης}
|
||||
\InstructorMail{rovithak@auth.gr}
|
||||
|
||||
\CoInstructorName{Θ. Αφορόζη}
|
||||
\CoInstructorMail{taforozi@ece.auth.gr}
|
||||
|
||||
\CurrentDate{\today}
|
||||
|
||||
|
||||
\usepackage{enumitem}
|
||||
\usepackage{tabularx}
|
||||
\usepackage{array}
|
||||
\usepackage{amssymb}
|
||||
\usepackage{amsfonts}
|
||||
\usepackage{amsmath}
|
||||
\usepackage{commath}
|
||||
|
||||
\usepackage{float}
|
||||
\usepackage[labelformat=empty]{subcaption}
|
||||
|
||||
\begin{document}
|
||||
|
||||
\InsertTitle
|
||||
|
||||
%\tableofcontents
|
||||
|
||||
\section{Εισαγωγή}
|
||||
Η παρούσα εργασία αφορά το πρόβλημα της ελαχιστοποίησης μιας δοσμένης συνάρτησης πολλών μεταβλητών $f: \mathbb{R}^n \rightarrow \mathbb{R}$ με περιορισμούς χρησιμοποιώντας τη μέθοδο μέγιστης καθόδου με προβολή.
|
||||
Η μέθοδος αυτή θα εκτελεστεί σε αντιπαραβολή με την αντίστοιχη μέθοδο χωρίς περιορισμούς από την προηγούμενη εργασία.
|
||||
Για το λόγο αυτό χρησιμοποιούμε των κώδικα της προηγούμενης εργασίας με κάποιες τροποποιήσεις, όπως θα δούμε και παρακάτω.
|
||||
|
||||
\subsection{Παραδοτέα}
|
||||
Τα παραδοτέα της εργασίας αποτελούνται από:
|
||||
\begin{itemize}
|
||||
\item Την παρούσα αναφορά.
|
||||
\item Τον κατάλογο \textbf{scripts/}, που περιέχει τον κώδικα της MATLAB.
|
||||
\item Το \href{https://git.hoo2.net/hoo2/OptimizationTechniques/src/branch/master/Work%203}{σύνδεσμο} με το αποθετήριο που περιέχει όλο το project με τον κώδικα της MATLAB, της αναφοράς και τα παραδοτέα.
|
||||
\end{itemize}
|
||||
|
||||
\subsection{Προγραμματιστική προσέγγιση}
|
||||
Για τον προγραμματισμό και εκτέλεση των μεθόδων της παρούσας εργασίας έγινε χρήση της MATLAB.
|
||||
Στον κατάλογο \textbf{scripts}, περιέχονται όλες οι μέθοδοι και οι τεχνικές υπολογισμού βημάτων με τη μορφή συναρτήσεων καθώς και scripts που τις καλούν.
|
||||
Για κάθε ένα θέμα της εργασίας, υπάρχει το αντίστοιχο script που περιέχει τους υπολογισμούς, τις κλήσεις των μεθόδων και τη δημιουργία των διαγραμμάτων.
|
||||
Για το πρώτο θέμα το αρχείο Script\_1\_SteepDesc.m για το δεύτερο το Script\_2\_SteepDesc\_Proj.m και ούτω καθεξής.
|
||||
Η μέθοδος μέγιστης καθόδου (αρχείο: \textbf{method\_SteepDesc.m}) είναι η ίδια με αυτή της προηγούμενης εργασίας με τη μόνη διαφορά ότι τροποποιήθηκε ώστε η αντικειμενική συνάρτηση να δέχεται διάνυσμα ώς όρισμα και όχι δύο διαφορετικές μεταβλητές $x, y$.
|
||||
Αυτό ακολουθήθηκε και για την έκδοση με προβολή και φυσικά είχε αντίκτυπο και στις υπόλοιπες συναρτήσεις, όπως η κλίση ή ο Εσσιανός.
|
||||
Στην παρούσα εργασία η υλοποίηση του κώδικα ακολουθεί την προσέγγιση των προηγούμενων εργασιών.
|
||||
Πιο συγκεκριμένα.
|
||||
|
||||
\subsection{Μέθοδοι επιλογής βήματος}
|
||||
Εφόσον στην προηγούμενη εργασία οι μέθοδοι επιλογής βήματος ήταν ανεξάρτητες από την μέθοδο υπολογισμού του ελάχιστου και εφόσον χρησιμοποιούμε τον ίδιο κώδικά και στην παρούσα εργασία, αυτός ο τρόπος σχεδίασης παρέμεινε.
|
||||
Ουσιαστικά για κάθε ένα τρόπο υπολογισμού του $\gamma_k$, υπάρχει αντίστοιχη συνάρτηση, με κοινό interface.
|
||||
Αυτό έχει τη μορφή: \textit{\textbf{gamma\_<method>(f, grad\_f, dk, xk)}}, όπου το \textbf{f} είναι η αντικειμενική συνάρτηση, \textbf{grad\_f} η συνάρτηση κλίσης της, \textbf{dk} η τιμή της συνάρτησης κλίσης στο xk και \textbf{xk} το σημείο ενδιαφέροντος.
|
||||
Έτσι οι μέθοδοι αντιγράφηκαν και εδώ για ολότητα, ακόμα και αν για την παρούσα εργασία χρησιμοποιείται μόνο το σταθερό βήμα $\gamma_k$.
|
||||
|
||||
%\subsection{Symbolic expression functions}
|
||||
%Μία ακόμη προγραμματιστική τεχνική που ακολουθήθηκε είναι η χρήση \textbf{symbolic expression} για την αναπαράσταση της αντικειμενικής συνάρτησης.
|
||||
%Η εξήγηση υπάρχει και στις προηγούμενες εργασίες αλλά την παραθέτουμε εδώ για ολότητα.
|
||||
%Ο λόγος που επιλέχθηκε είναι η \textbf{δυνατότητα εξαγωγής ενός symbolic expression που αναπαριστά την κλίση $\nabla f$ και τον Εσσιανό $\nabla^2f$ μιας συνάρτησης} από την MATLAB, κάνοντας χρήση των εντολών \textit{gradient()} και \textit{hessian()}.
|
||||
%Αν αντίθετα χρησιμοποιούσαμε απλές συναρτήσεις, πολυώνυμα ή lambdas για την αναπαράσταση των αντικειμενικών συναρτήσεων, τότε για τον υπολογισμό της κλίσης και του Εσσιανού θα έπρεπε:
|
||||
%\begin{itemize}
|
||||
% \item Είτε να υπολογίζαμε αριθμητικά τις παραγώγους gradient και hessian μέσα στις μεθόδους, κάτι που θα εισήγαγε \textit{\textbf{αχρείαστο αριθμητικό σφάλμα}}.
|
||||
% \item Είτε να κάναμε χρήση δύο επιπλέων συναρτήσεων (ή πολυωνύμων) για την αναπαράσταση τους, κάτι που ουσιαστικά θα δημιουργούσε \textit{\textbf{πλεονασμό πληροφορίας εισόδου}} και άρα μεγαλύτερη πιθανότητα να κάνουμε λάθος.
|
||||
%\end{itemize}
|
||||
%Η αναπαράσταση όμως με χρήση symbolic expression είναι πιο “βαριά” όταν χρειάζεται να υπολογίσουμε την τιμή μιας συνάρτησης σε κάποιο σημείο (subs(expr, number)).
|
||||
%Αυτό είναι κάτι που χρειάζεται εκτενώς στον κώδικά μας.
|
||||
%Για το λόγο αυτό, ενώ η συνάρτηση δίνεται ως symbolic expression, μέσω αυτής υπολογίζονται αυτόματα η κλίση, ο Εσσιανός αλλά και οι “κανονικές” συναρτήσεις MATLAB που τις υλοποιούν.
|
||||
%Έτσι έχουμε την ακριβή αναπαράσταση της κλίσης και του Εσσιανού ως συναρτήσεις χωρίς να πληρώνουμε το κόστος της subs().
|
||||
|
||||
\section{Απεικόνιση της συνάρτησης}
|
||||
Η συνάρτηση με την οποία ασχολούμαστε στην παρούσα εργασία είναι η:
|
||||
\boldmath
|
||||
\begin{equation}
|
||||
f: \mathbb{R}^2 \rightarrow \mathbb{R}, f(x) = \frac{1}{3}{x_1}^2 + 3{x_2}^2, x = \begin{bmatrix} x_1 \\ x_2 \end{bmatrix}
|
||||
\end{equation}
|
||||
\label{eq:ObjectiveFunction}
|
||||
Όπου $ -10 \leq x_1 \leq 5 $ και $ -8 \leq x_2 \leq 12 $ \\
|
||||
\unboldmath
|
||||
Στο παρακάτω σχήμα \ref{fig:plot3dFunction} φαίνεται η τρισδιάστατη απεικόνιση της συνάρτησης.
|
||||
\InsertFigure{!h}{0.8}{fig:plot3dFunction}{../scripts/figures/Plot_Function.png}{Γραφική παράσταση της f}
|
||||
|
||||
Από το σχήμα μπορούμε πολύ εύκολα να διακρίνουμε ότι η συνάρτηση είναι κυρτή στο σύνολο των περιορισμών της εκφώνησης $ -10 \leq x_1 \leq 5 $ και $ -8 \leq x_2 \leq 12 $.
|
||||
Για να πάρουμε μια καλύτερη αίσθηση για το που βρίσκεται το τοπικό ελάχιστο της $f$, παρακάτω παραθέτουμε ένα γράφημα με τις ισοβαρείς καμπύλες της $f$.
|
||||
\InsertFigure{H}{0.8}{fig:plotContour}{../scripts/figures/Plot_Contour.png}{Ισοβαρείς της f}
|
||||
|
||||
Από το παραπάνω σχήμα \ref{fig:plotContour} φαίνονται και γραφικά οι μικρές κλίσης που παρουσιάζει η συνάρτηση κοντά στο ελάχιστο σημείο $(0,0)$.
|
||||
Τα διαγράμματα για τη μέθοδο δημιουργούνται εκτελώντας το αρχείο \textbf{Script\_0\_Plots.m}
|
||||
|
||||
\section{Θέμα 1 - Μέθοδος Μέγιστης Καθόδου χωρίς περιορισμούς}
|
||||
Εφαρμόζοντας την μέθοδο μέγιστης καθόδου από την προηγούμενη εργασία, με ακρίβεια $\epsilon = 0.001$, για τα βήματα $\gamma_k$ της εκφώνησης, παρατηρούμε ότι η μέθοδος συγκλίνει στο ελάχιστο για μικρά $\gamma_k$ ενώ \textbf{αποκλίνει για μεγάλα} \boldmath$\gamma_k > 0.33$\unboldmath.
|
||||
Από τις δοκιμές φαίνεται ότι το σημείο εκκίνησης δεν παίζει ρόλο και για αυτό επιλέξαμε να παραθέσουμε τα ευρήματά μας από το σημείο $(5,-5)$, για αντιπαραβολή με το επόμενο βήμα της εκφώνησης.
|
||||
\InsertFigure{H}{0.6}{fig:StDes_Iter_o_gamma_2}{../scripts/figures/StDes_Iter_o_gamma_1.png}{Αριθμός επαναλήψεων για διαφορετικές τιμές $\gamma_k$ [Μέγιστη Κάθοδος].}
|
||||
Επίσης παρατηρούμε ότι για μικρό \boldmath$\gamma_k = 0.1$ η σύγκλιση είναι ομαλή, ενώ για μεγάλο $\gamma_k = 0.3$ \unboldmath παρουσιάζει ταλάντωση κατά την σύγκλιση.
|
||||
Παρακάτω στο σχήμα \ref{fig:StDes_gamma} παραθέτουμε την πορεία σύγκλισης και απόκλισης για τις διαφορετικές τιμές του $\gamma_k$.
|
||||
|
||||
\begin{figure}[ht]
|
||||
\centering
|
||||
% First row
|
||||
\begin{subfigure}{0.45\textwidth}
|
||||
\centering
|
||||
\includegraphics[width=\linewidth]{../scripts/figures/StDes_gamma_0.1.png}
|
||||
\caption{$\gamma_k = 0.1$}
|
||||
\label{fig:StDes_gamma_0.1}
|
||||
\end{subfigure}
|
||||
\hfill
|
||||
\begin{subfigure}{0.45\textwidth}
|
||||
\centering
|
||||
\includegraphics[width=\linewidth]{../scripts/figures/StDes_gamma_0.3.png}
|
||||
\caption{$\gamma_k = 0.3$}
|
||||
\label{fig:StDes_gamma_0.3}
|
||||
\end{subfigure}
|
||||
|
||||
% Second row
|
||||
\vspace{1em}
|
||||
\begin{subfigure}{0.45\textwidth}
|
||||
\centering
|
||||
\includegraphics[width=\linewidth]{../scripts/figures/StDes_gamma_3.png}
|
||||
\caption{$\gamma_k = 3$}
|
||||
\label{fig:StDes_gamma_3}
|
||||
\end{subfigure}
|
||||
\hfill
|
||||
\begin{subfigure}{0.45\textwidth}
|
||||
\centering
|
||||
\includegraphics[width=\linewidth]{../scripts/figures/StDes_gamma_5.png}
|
||||
\caption{$\gamma_k = 5$}
|
||||
\label{fig:StDes_gamma_5}
|
||||
\end{subfigure}
|
||||
|
||||
\caption{Σύγκριση της μεθόδου Steepest descent για διαφορετικά $\gamma_k$}
|
||||
\label{fig:StDes_gamma}
|
||||
\end{figure}
|
||||
|
||||
\subsection{Μαθηματική ανάλυση}
|
||||
Τα παραπάνω αποτελέσματα επιβεβαιώνονται και θεωρητικά.
|
||||
Πιο συγκεκριμένα για τη σύγκλιση της μεθόδου μέγιστης καθόδου όπου το κάθε σημείο υπολογίζεται από την σχέση \boldmath$x_{k+1} = x_k - \gamma_k \nabla f(x_k)$, \unboldmath πρέπει να ισχύουν:
|
||||
\begin{enumerate}
|
||||
\item H $f$ να είναι κυρτή.
|
||||
\item Η $f$ να είναι συνεχής και διαφορίσιμη και η κλίση της υπολογίσιμη.
|
||||
\item Για το βήμα υπολογισμού να ισχύει η σχέση:
|
||||
\begin{equation} 0 < \gamma_k < \frac{2}{L} \label{eq:gammaLimmit} \end{equation}
|
||||
Όπου $L$ το άνω φράγμα της Lipschitz για την κλίση $\nabla f(x)$ (αν είναι γνωστή), η οποία είναι η μέγιστη ιδιοτιμή του Εσσιανού και δίνεται από τη σχέση:
|
||||
\begin{equation}
|
||||
L = \max_{x} \{\lambda_{max} (H(x))\}
|
||||
\label{eq:Lipschitz}
|
||||
\end{equation}
|
||||
\end{enumerate}
|
||||
\par
|
||||
Έτσι για τη δική μας περίπτωση η κλίση της $f(x)$ είναι:
|
||||
\[
|
||||
\nabla f(x) = \begin{bmatrix} \frac{\partial f}{\partial x1} \\ \frac{\partial f}{\partial x2} \end{bmatrix} = \begin{bmatrix} \frac{2}{3}x_1 \\ 6x_2 \end{bmatrix}
|
||||
\]
|
||||
Ο Εσσιανός πίνακας $H(x)$ της $f(x)$ είναι:
|
||||
\[
|
||||
H(x) =
|
||||
\begin{bmatrix}
|
||||
\frac{\partial^2 f}{\partial x_1^2} & \frac{\partial^2 f}{\partial x_1 \partial x_2} \\
|
||||
\frac{\partial^2 f}{\partial x_2 \partial x_1} & \frac{\partial^2 f}{\partial x_2^2}
|
||||
\end{bmatrix}
|
||||
=
|
||||
\begin{bmatrix}
|
||||
\frac{2}{3} & 0 \\ 0 & 6
|
||||
\end{bmatrix}
|
||||
\]
|
||||
Και εφόσον είναι διαγώνιος οι ιδιοτιμές του είναι οι τιμές της διαγωνίου. Δηλαδή:
|
||||
\[
|
||||
\det(H - \lambda I) = 0 \Leftrightarrow
|
||||
\det(\begin{bmatrix}\frac{2}{3} - \lambda & 0 \\ 0 & 6 - \lambda \end{bmatrix}) = 0 \Leftrightarrow
|
||||
\left(\frac{2}{3} - \lambda\right)(6 - \lambda) = 0 \Leftrightarrow
|
||||
\lambda_{\min} = \frac{2}{3}, \quad \lambda_{\max} = 6.
|
||||
\]
|
||||
Έτσι από τις εξισώσεις (\ref{eq:gammaLimmit}) και (\ref{eq:Lipschitz}) προκύπτει τελικά ότι $\frac{2}{L} = \frac{2}{\lambda_{max}} = \frac{1}{3}$ και επομένως για σύγκλιση της μεθόδου πρέπει να ισχύει:
|
||||
\boldmath
|
||||
\begin{equation}
|
||||
0 < \gamma_k < \frac{1}{3}
|
||||
\label{eq:gammaConvergence}
|
||||
\end{equation}
|
||||
\unboldmath
|
||||
Βλέπουμε ότι από την ανάλυσή μας, η τιμή που βρήκαμε εμπειρικά για την απόκλιση από την εκτέλεση του αλγορίθμου για το $\gamma > 0.33$ επιβεβαιώνεται.
|
||||
\par
|
||||
Επίσης το γεγονός ότι ο Εσσιανός είναι διαγώνιος, μας δίνει επίσης πληροφορία ότι:
|
||||
\boldmath
|
||||
\begin{itemize}
|
||||
\item Η ιδιοτιμή $\lambda_1 = \frac{2}{3}$ εφαρμοσμένη στη σχέση (\ref{eq:gammaLimmit}) μας δίνει: \begin{equation} 0<\gamma_k<3 \label{eq:gammaLimmit_x1} \end{equation}
|
||||
που είναι η τιμή για την οποία η μέθοδος συγκλίνει στη διάσταση $x_1$
|
||||
\item Η ιδιοτιμή $\lambda_2 = 6$ εφαρμοσμένη στη σχέση (\ref{eq:gammaLimmit}) μας δίνει:
|
||||
\begin{equation} 0<\gamma_k<\frac{1}{3}\label{eq:gammaLimmit_x2} \end{equation}
|
||||
που είναι η τιμή για την οποία η μέθοδος συγκλίνει στη διάσταση $x_2$
|
||||
\end{itemize}
|
||||
\unboldmath
|
||||
|
||||
\par
|
||||
\underline{Εναλλακτικά} \\[0.5ex]
|
||||
Αν θέλαμε να βρούμε το κριτήριο σύγκλισης για το βήμα $\gamma_k$ ξεχωριστά για την κάθε διάσταση θα μπορούσαμε να θεωρήσουμε ισοδύναμα με την παραπάνω ανάλυση, πως για να συγκλίνει η μέθοδος θα πρέπει να ισχύει:
|
||||
\[
|
||||
\abs{\frac{x_{k+1}}{x_k}} < 1
|
||||
\]
|
||||
Από την παραπάνω εξίσωση προκύπτει πάλι ένα σύστημα εξισώσεων για την κάθε διάσταση:
|
||||
\[
|
||||
\left\{
|
||||
\begin{aligned}
|
||||
\abs{\dfrac{x_{1,k} - \frac{2}{3}\gamma_k x_{1,k}}{x_{1,k}}} < 1 \\
|
||||
\abs{\dfrac{x_{2,k} - 6\gamma_k x_{2,k}}{x_{2,k}}} < 1
|
||||
\end{aligned}
|
||||
\right.
|
||||
\Leftrightarrow
|
||||
\left\{
|
||||
\begin{aligned}
|
||||
\abs{1 - \frac{2}{3}\gamma_k} < 1 \\
|
||||
\abs{1 - 6\gamma_k} < 1
|
||||
\end{aligned}
|
||||
\right.
|
||||
\Leftrightarrow
|
||||
\left\{
|
||||
\begin{aligned}
|
||||
0 < \gamma_k < 3 \\
|
||||
0 < \gamma_k < \frac{1}{3}
|
||||
\end{aligned}
|
||||
\right.
|
||||
\]
|
||||
Που επιβεβαιώνει την προηγούμενη ανάλυσή μας για τη συνολική σύγκλιση, αλλά και για την κάθε διάσταση ξεχωριστά.
|
||||
|
||||
\section{Μέθοδος Μέγιστης Καθόδου με προβολή}
|
||||
Πριν περάσουμε στις υπόλοιπες απαιτήσεις της εργασίας θα θέλαμε να παραθέσουμε κάποιες πληροφορίες για την υλοποίηση της μεθόδου μέγιστης καθόδου με προβολή (αρχείο: \textbf{method\_SteepDesc\_Proj.m}).
|
||||
Η συνάρτηση αυτή δέχεται ως είσοδο την αντικειμενική συνάρτηση και την συνάρτηση κλίσης καθώς και το σημείο εκκίνησης $x_k$ και το βήμα $s_k$.
|
||||
Με τη βοήθεια της συνάρτησης \textit{ProjectionPoint()} \textbf{παίρνει πρώτα την προβολή} του $x_k$ στο διάστημα των περιορισμών αν αυτό χρειάζεται \textbf{και έπειτα εφαρμόζει τον αλγόριθμο}.
|
||||
Αυτό σημαίνει ότι μπορεί να χρησιμοποιηθεί και για σημεία εκκίνησης εκτός του συνόλου των περιορισμών.
|
||||
Ο αλγόριθμος είναι παρόμοιος με αυτόν της προηγούμενης εργασίας με τη διαφορά ότι η διεύθυνση $d_k$ επιλέγεται από τη σχέση:
|
||||
\[
|
||||
d_k = Pr_X\{ x_k - s_k \nabla f(x_k)\} - x_k = \bar{x_k} - x_k
|
||||
\]
|
||||
και τελικά έχουμε:
|
||||
\[
|
||||
x_{k+1} = x_k + \gamma_k d_k
|
||||
\]
|
||||
Δηλαδή εφαρμόζουμε πρώτα τη μέθοδο μέγιστης καθόδου με βήμα $s_k$ στην κατεύθυνση $-\nabla f$ και έπειτα προβάλουμε το σημείο στο σύνολο $X$ και χρησιμοποιούμε αυτό ως διεύθυνση με βήμα $\gamma_k$.
|
||||
\par
|
||||
Αξίζει να παρατηρήσουμε πως \textbf{αν το σημείο $x_k$ είναι εντός του συνόλου περιορισμών} τότε η προβολή του σημείου στο σύνολο, είναι το ίδιο το σημείο και έτσι:
|
||||
\[
|
||||
\begin{aligned}
|
||||
x_{k+1} & = x_k + \gamma_k \left( Pr_X\{ x_k - s_k \nabla f(x_k)\} - x_k \right) \\
|
||||
& = x_k + \gamma_k \left( ( x_k - s_k \nabla f(x_k) ) - x_k \right) \\
|
||||
& = x_k - \gamma_k s_k \nabla f(x_k) \\
|
||||
& = x_k - \gamma_k' \nabla f(x_k)
|
||||
\end{aligned}
|
||||
\]
|
||||
Βλέπουμε λοιπόν ότι για σημεία εσωτερικά του Χ, έχουμε ουσιαστικά τη μέθοδο της μέγιστης καθόδου με βήμα $\gamma_k' = \gamma_k s_k$.
|
||||
Εφόσον όμως για μη εφικτά σημεία, πριν εφαρμόσουμε τη μέθοδο, τα προβάλουμε στο σύνολο των περιορισμών Χ, βλέπουμε πως η \textbf{μέθοδος τελικά προσεγγίζει τη μέγιστη κάθοδο με βήμα $\gamma_k'$, που είναι και το βήμα για το οποίο ελέγχουμε αν ισχύουν τα κριτήρια σύγκλισης}.
|
||||
|
||||
\section{Θέμα 2 - Σημείο (5, -5), $s_k = 5, \gamma_k = 0.5$}
|
||||
Εφαρμόζοντας τη μέθοδο για ακρίβεια $\epsilon = 0.01$, $s_k = 5$ και $\gamma_k = 0.5$ έχουμε $\gamma_k' = 2.5 > \frac{1}{3}$, άρα το κριτήριο σύγκλισης δεν πληρείται.
|
||||
\InsertFigure{H}{0.8}{fig:StDesProj_sk_5_gamma_0.5}{../scripts/figures/StDesProj_sk_5_gamma_0.5.png}{Μέθοδος μέγιστης καθόδου με προβολή για $s_k = 5, \gamma_k = 0.5$.}
|
||||
|
||||
Παρατηρούμε πως ενώ η μέθοδος ταλαντώνει και δεν συγκλίνει στο ελάχιστο \textbf{όπως και η αντίστοιχη εκτέλεση της μέγιστης καθόδου χωρίς περιορισμούς με το ίδιο} \boldmath$\gamma_k$\unboldmath.
|
||||
Παρόλα αυτά όμως, η ταλάντωση λαμβάνει χώρα \textbf{εντός του συνόλου των περιορισμών} της εκφώνησης.
|
||||
|
||||
\section{Θέμα 3 - Σημείο (-5, 10), $s_k = 15, \gamma_k = 0.1$}
|
||||
Εφαρμόζοντας τη μέθοδο για ακρίβεια $\epsilon = 0.01$, $s_k = 15$ και $\gamma_k = 0.1$ έχουμε $\gamma_k' = 1.5 > \frac{1}{3}$, άρα το κριτήριο σύγκλισης δεν πληρείται και πάλι.
|
||||
\InsertFigure{H}{0.8}{fig:StDesProj_sk_15_gamma_0.1}{../scripts/figures/StDesProj_sk_15_gamma_0.1.png}{Μέθοδος μέγιστης καθόδου με προβολή για $s_k = 15, \gamma_k = 0.1$.}
|
||||
|
||||
Και εδώ παρατηρούμε πως ενώ το $\gamma_k$ έχει επιλεγεί θεωρητικά στο εύρος που οδηγεί σε σύγκλιση, το αντίστοιχο βήμα $s_k$ είναι πολύ μεγάλο, με αποτέλεσμα το γινόμενό τους $\gamma_k' = 15 * 0.1 = 1.5$ συνολικά να μην πληροί το κριτήριο $\gamma_k' < \frac{1}{3}$ και η μέθοδος να ταλαντώνει και πάλι.
|
||||
Αυτή τη φορά μόνο στον άξονα $x_2$.
|
||||
Αυτό βέβαια εξηγείται από την παραπάνω ανάλυσή καθώς βλέπουμε πως ενώ το γινόμενο 1.5 δεν πληροί τις προϋποθέσεις για σύγκλιση, πληροί όμως τις προϋποθέσεις της εξίσωσης (\ref{eq:gammaLimmit_x1}) με αποτέλεσμα να ταλαντώνει μόνο στον άξονα $x_2$
|
||||
Αυτό φυσικά είναι αληθές και για την προηγούμενη περίπτωση όπου το $\gamma_k' = 5 * 0.5 = 2.5$ και όπου πάλι η σύγκλιση ήταν μερική, μόνο για την διάσταση $x_1$.
|
||||
|
||||
\section{Θέμα 4 - Σημείο (8, -10), $s_k = 0.1, \gamma_k = 0.2$ - Θέμα 4}
|
||||
Αρχικά παρατηρούμε πως το σημείο \textbf{δεν είναι εφικτό}, καθώς είναι εκτός του συνόλου των περιορισμών της εκφώνησης.
|
||||
Αυτό βέβαια δεν μας αποτρέπει από την εφαρμογή της μεθόδου, καθώς αρχικά μπορούμε να προβάλουμε το σημείο στο σύνολο και να εφαρμόσουμε τη μέθοδο έπειτα.
|
||||
Ακόμα, αυτή τη φορά οι τιμές των βημάτων $s_k, \gamma_k$, έχουν επιλεγεί μέσα στο εύρος για το οποίο έχουμε σύγκλιση, καθώς: $\gamma_k' = 0.1 * 0.2 = 0.02 < \frac{1}{3}$, επομένως αναμένουμε η μέθοδος να συγκλίνει στο ελάχιστο.
|
||||
Εφαρμόζοντας τη μέθοδο για ακρίβεια $\epsilon = 0.01$ έχουμε:
|
||||
\InsertFigure{H}{0.8}{fig:StDesProj_sk_0.1_gamma_0.2}{../scripts/figures/StDesProj_sk_0.1_gamma_0.2.png}{Μέθοδος μέγιστης καθόδου με προβολή για $s_k = 0.1, \gamma_k = 0.2$.}
|
||||
Όπου βλέπουμε πως η μέθοδος συγκλίνει επιβεβαιώνοντας την παραπάνω ανάλυση.
|
||||
|
||||
\section{Συμπεράσματα}
|
||||
|
||||
Η μέθοδος μέγιστης καθόδου, με και χωρίς προβολή, παρουσιάζει τις αναμενόμενες συμπεριφορές σύγκλισης ανάλογα με την επιλογή του βήματος $\gamma_k$.
|
||||
Για την περίπτωση χωρίς προβολή, επιβεβαιώθηκε θεωρητικά και εμπειρικά ότι η μέθοδος συγκλίνει μόνο όταν το $\gamma_k$ βρίσκεται εντός του εύρους που ορίζεται από την ανάλυση της Lipschitz σταθεράς, δηλαδή $0 < \gamma_k < \frac{1}{3}$.
|
||||
Επιπλέον, παρατηρήθηκε πως για μεγαλύτερες τιμές$\gamma_k$ η μέθοδος αποκλίνει, όπως προβλέπεται από την μαθηματική ανάλυση.
|
||||
|
||||
Αντίστοιχα, για τη μέθοδο με προβολή, αποδεικνύεται πως το κριτήριο σύγκλισης εξαρτάται από τον συνδυασμό του $s_k$ (βήμα αντίθετο στην κατεύθυνση της κλίσης) και του $\gamma_k$ (βήμα κατά την προβολή).
|
||||
Η ανάλυση επιβεβαίωσε ότι η μέθοδος προσεγγίζει τη μέθοδο χωρίς προβολή για εφικτά σημεία, ενώ για μη εφικτά σημεία εξασφαλίζει ότι οι διαδοχικές προσεγγίσεις παραμένουν εντός των περιορισμών.
|
||||
|
||||
Τέλος, μέσα από τα παραδείγματα, αναδείχθηκε η σημασία της σωστής επιλογής των παραμέτρων $\gamma_k$ και $s_k$ για τη συνολική σταθερότητα και σύγκλιση της μεθόδου, καθώς και ο τρόπος με τον οποίο η ανάλυση των ιδιοτιμών του Εσσιανού επιτρέπει την καλύτερη κατανόηση της συμπεριφοράς της μεθόδου.
|
||||
|
||||
\end{document}
|
||||
@@ -0,0 +1,19 @@
|
||||
function [flag] = BelongsTo(x, SetLimmits)
|
||||
%Checks if the x vector belongs to Set
|
||||
%
|
||||
% x: A vector to project
|
||||
% SetLimmits: The set to project. Each line/dimension off the set has to contain the limits
|
||||
% of the set to that particular dimension
|
||||
% flag: True if belongs
|
||||
%
|
||||
flag = true; % Have faith
|
||||
|
||||
for i = 1:size(SetLimmits, 2)
|
||||
if x(i) < SetLimmits(i,1)
|
||||
flag = false;
|
||||
elseif x(i) > SetLimmits(i,2)
|
||||
flag = false;
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
@@ -0,0 +1,36 @@
|
||||
% Given environment
|
||||
|
||||
clear;
|
||||
% Setup the function under test
|
||||
syms x [2 1] real;
|
||||
fexpr = (1/3)*x(1)^2 +3*x(2)^2;
|
||||
title_fun = "$f(x) = \frac{1}{3}{x_1}^2 + 3{x_2}^2$";
|
||||
|
||||
XSetLimmits = [-10, 5 ; -8, 12];
|
||||
|
||||
% Calculate the gradient and Hessian
|
||||
grad_fexpr = gradient(fexpr, x); % Gradient of f
|
||||
hessian_fexpr = hessian(fexpr, x); % Hessian of f
|
||||
|
||||
% Convert symbolic expressions to MATLAB functions
|
||||
fun = matlabFunction(fexpr, 'Vars', {x}); % Function
|
||||
grad_fun = matlabFunction(grad_fexpr, 'Vars', {x}); % Gradient
|
||||
hessian_fun = matlabFunction(hessian_fexpr, 'Vars', {x}); % Hessian
|
||||
|
||||
% Minimum reference
|
||||
[Xmin, Fmin] = fminsearch(fun, [-1, -1]');
|
||||
Xmin = round(Xmin, 3);
|
||||
Fmin = round(Fmin, 3);
|
||||
|
||||
% 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])
|
||||
|
||||
%fixed step size globals
|
||||
global gamma_fixed_step
|
||||
|
||||
global image_width,
|
||||
global image_height;
|
||||
|
||||
image_width = 960;
|
||||
image_height = 640;
|
||||
@@ -0,0 +1,18 @@
|
||||
function [x_p] = ProjectionPoint(x, SetLimmits)
|
||||
%Project the x vector to Set, returns a point of Set close(st) to x
|
||||
%
|
||||
% x: A vector to project
|
||||
% SetLimmits: The set to project. Each line/dimension off the set has to contain the limits
|
||||
% of the set to that particular dimension
|
||||
%x_p: The projected point
|
||||
%
|
||||
x_p = x;
|
||||
for i = 1:size(SetLimmits, 2)
|
||||
if x(i) < SetLimmits(i,1)
|
||||
x_p(i) = SetLimmits(i,1);
|
||||
elseif x(i) > SetLimmits(i,2)
|
||||
x_p(i) = SetLimmits(i,2);
|
||||
end
|
||||
end
|
||||
end
|
||||
|
||||
@@ -0,0 +1,14 @@
|
||||
% Define environment (functions, gradients etc...)
|
||||
GivenEnv
|
||||
|
||||
%
|
||||
% We plot the function in the domain of x,y in [-3, 3].
|
||||
% We also plot the contour in order to get a sense of the min and maximum
|
||||
% points in the x-y plane
|
||||
%
|
||||
|
||||
% 3d plot the function
|
||||
plot3dFun(fun, XSetLimmits(1, :), XSetLimmits(2, :), 100, title_fun, "figures/Plot_Function.png");
|
||||
|
||||
% Plot isobaric lines
|
||||
plotContour(fun, XSetLimmits(1, :), XSetLimmits(2, :), 100, title_fun, "figures/Plot_Contour.png");
|
||||
@@ -0,0 +1,40 @@
|
||||
% Define environment (functions, gradients etc...)
|
||||
GivenEnv
|
||||
|
||||
% Define parameters
|
||||
max_iter = 1000; % Maximum iterations
|
||||
tol = 0.001; % Tolerance
|
||||
|
||||
% Point x0 = (1, 1)
|
||||
% =========================================================================
|
||||
point = 1;
|
||||
x0 = [5, -5]';
|
||||
point_str = "[" + x0(1) + ", " + x0(2) + "]";
|
||||
f = fun(x0);
|
||||
gf = grad_fun(x0);
|
||||
hf = hessian_fun(x0);
|
||||
fprintf('Initial point (%d, %d), f = %f, grad = [%f;%f], hessian = [%f %f ; %f %f]=> Method applicable\n', x0, f, gf, hf);
|
||||
|
||||
|
||||
k = zeros(100, 1);
|
||||
j = 1;
|
||||
n = linspace(0.1, 1.5, 100);
|
||||
for g = n
|
||||
gamma_fixed_step = g;
|
||||
[~, ~, k(j)] = method_SteepDesc(fun, grad_fun, x0, tol, max_iter, 'fixed');
|
||||
j = j + 1;
|
||||
end
|
||||
plotItersOverGamma(n, k, "Iteration for different $\gamma$ values", "figures/StDes_Iter_o_gamma_" + point + ".png");
|
||||
|
||||
for g=[0.1, 0.3, 3, 5]
|
||||
gamma_fixed_step = g;
|
||||
|
||||
[x_fixed, f_fixed, kk] = method_SteepDesc(fun, grad_fun, x0, tol, max_iter, 'fixed');
|
||||
fprintf('Fixed step g=%f: Initial point (%f, %f), steps:%d, Final (x1,x2)=(%f, %f), f(x1,x2)=%f\n', g, x0, kk, x_fixed(:, end), f_fixed(end));
|
||||
if g <= 1
|
||||
plotPointsOverContour(x_fixed, fun, XSetLimmits(1, :), XSetLimmits(2, :), 100, point_str + ": Steepest descent $\gamma$ = " + gamma_fixed_step, "figures/StDes_gamma_" + g + ".png");
|
||||
else
|
||||
plotPointsOverContour(x_fixed, fun, [-500, 500], [-500, 500], 100, point_str + ": Steepest descent $\gamma$ = " + gamma_fixed_step, "figures/StDes_gamma_" + g + ".png");
|
||||
end
|
||||
end
|
||||
|
||||
@@ -0,0 +1,32 @@
|
||||
% Define environment (functions, gradients etc...)
|
||||
GivenEnv
|
||||
|
||||
% Settings
|
||||
max_iter = 1000; % Maximum iterations
|
||||
|
||||
% Define parameters
|
||||
% =========================================================================
|
||||
x0 = [5, -5]';
|
||||
gamma = 0.5;
|
||||
sk_step = 5;
|
||||
tol = 0.01; % Tolerance
|
||||
|
||||
|
||||
|
||||
point_str = "[" + x0(1) + ", " + x0(2) + "]";
|
||||
f = fun(x0);
|
||||
gf = grad_fun(x0);
|
||||
hf = hessian_fun(x0);
|
||||
fprintf('Initial point (%d, %d), f = %f, grad = [%f;%f], hessian = [%f %f ; %f %f]=> Method applicable\n', x0, f, gf, hf);
|
||||
|
||||
|
||||
gamma_fixed_step = gamma;
|
||||
[x_fixed, f_fixed, kk] = method_SteepDesc_Proj(fun, grad_fun, x0, sk_step, XSetLimmits, tol, max_iter, 'fixed');
|
||||
fprintf('Fixed step g=%f: Initial point (%f, %f), steps:%d, Final (x1,x2)=(%f, %f), f(x1,x2)=%f\n', gamma_fixed_step, x0, kk, x_fixed(:, end), f_fixed(end));
|
||||
plotPointsOverContour(x_fixed, fun, XSetLimmits(1, :), XSetLimmits(2, :), 100, point_str + ": Steepest descent proj $s_k$ = " + sk_step + ", $\gamma$ = " + gamma_fixed_step, "figures/StDesProj_sk_" + sk_step + "_gamma_" + gamma + ".png");
|
||||
|
||||
%[x_minimized, f_minimized, kk] = method_SteepDesc_Proj(fun, grad_fun, x0, sk_step, XSetLimmits, tol, max_iter, 'minimized');
|
||||
%fprintf('Minimized f: Initial point (%f, %f), steps:%d, Final (x1,x2)=(%f, %f), f(x1,x2)=%f\n', x0, kk, x_fixed(:, end), f_fixed(end));
|
||||
%plotPointsOverContour(x_fixed, fun, XSetLimmits(1, :), XSetLimmits(2, :), 100, point_str + ": Steepest descent minimized", "");
|
||||
|
||||
|
||||
@@ -0,0 +1,27 @@
|
||||
% Define environment (functions, gradients etc...)
|
||||
GivenEnv
|
||||
|
||||
% Settings
|
||||
max_iter = 1000; % Maximum iterations
|
||||
|
||||
% Define parameters
|
||||
% =========================================================================
|
||||
x0 = [-5, 10]';
|
||||
gamma = 0.1;
|
||||
sk_step = 15;
|
||||
tol = 0.01; % Tolerance
|
||||
|
||||
|
||||
|
||||
point_str = "[" + x0(1) + ", " + x0(2) + "]";
|
||||
f = fun(x0);
|
||||
gf = grad_fun(x0);
|
||||
hf = hessian_fun(x0);
|
||||
fprintf('Initial point (%d, %d), f = %f, grad = [%f;%f], hessian = [%f %f ; %f %f]=> Method applicable\n', x0, f, gf, hf);
|
||||
|
||||
|
||||
gamma_fixed_step = gamma;
|
||||
[x_fixed, f_fixed, kk] = method_SteepDesc_Proj(fun, grad_fun, x0, sk_step, XSetLimmits, tol, max_iter, 'fixed');
|
||||
fprintf('Fixed step g=%f: Initial point (%f, %f), steps:%d, Final (x1,x2)=(%f, %f), f(x1,x2)=%f\n', gamma_fixed_step, x0, kk, x_fixed(:, end), f_fixed(end));
|
||||
plotPointsOverContour(x_fixed, fun, XSetLimmits(1, :), XSetLimmits(2, :), 100, point_str + ": Steepest descent proj $s_k$ = " + sk_step + ", $\gamma$ = " + gamma_fixed_step, "figures/StDesProj_sk_" + sk_step + "_gamma_" + gamma + ".png");
|
||||
|
||||
@@ -0,0 +1,27 @@
|
||||
% Define environment (functions, gradients etc...)
|
||||
GivenEnv
|
||||
|
||||
% Settings
|
||||
max_iter = 1000; % Maximum iterations
|
||||
|
||||
% Define parameters
|
||||
% =========================================================================
|
||||
x0 = [8, -10]';
|
||||
gamma = 0.2;
|
||||
sk_step = 0.1;
|
||||
tol = 0.01; % Tolerance
|
||||
|
||||
|
||||
|
||||
point_str = "[" + x0(1) + ", " + x0(2) + "]";
|
||||
f = fun(x0);
|
||||
gf = grad_fun(x0);
|
||||
hf = hessian_fun(x0);
|
||||
fprintf('Initial point (%d, %d), f = %f, grad = [%f;%f], hessian = [%f %f ; %f %f]=> Method applicable\n', x0, f, gf, hf);
|
||||
|
||||
|
||||
gamma_fixed_step = gamma;
|
||||
[x_fixed, f_fixed, kk] = method_SteepDesc_Proj(fun, grad_fun, x0, sk_step, XSetLimmits, tol, max_iter, 'fixed');
|
||||
fprintf('Fixed step g=%f: Initial point (%f, %f), steps:%d, Final (x1,x2)=(%f, %f), f(x1,x2)=%f\n', gamma_fixed_step, x0, kk, x_fixed(:, end), f_fixed(end));
|
||||
plotPointsOverContour(x_fixed, fun, XSetLimmits(1, :), XSetLimmits(2, :), 100, point_str + ": Steepest descent proj $s_k$ = " + sk_step + ", $\gamma$ = " + gamma_fixed_step, "figures/StDesProj_sk_" + sk_step + "_gamma_" + gamma + ".png");
|
||||
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 91 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 1.1 MiB |
Binary file not shown.
|
After Width: | Height: | Size: 131 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 124 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 109 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 135 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 66 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 110 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 125 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 130 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 132 KiB |
@@ -0,0 +1,49 @@
|
||||
function [a, b, k, n] = fmin_bisection(fun, alpha, beta, epsilon, lambda)
|
||||
% Bisection method for finding the local minimum of a function.
|
||||
%
|
||||
% fun: The objective function
|
||||
% alpha: (number) The starting point of the interval in which we seek
|
||||
% for minimum
|
||||
% beta: (number) The ending point of the interval in which we seek
|
||||
% for minimum
|
||||
% epsilon: (number) The epsilon value (distance from midpoint)
|
||||
% lambda: (number) The lambda value (accuracy)
|
||||
%
|
||||
% return:
|
||||
% a: (vector) Starting points of the interval for each iteration
|
||||
% b: (vector) Ending points of the interval for each iteration
|
||||
% k: (number) The number of iterations
|
||||
% n: (number) The calls of objective function fun_expr
|
||||
%
|
||||
|
||||
|
||||
|
||||
% Error checking
|
||||
if alpha > beta || 2*epsilon >= lambda || lambda <= 0
|
||||
error ('Input criteria not met')
|
||||
end
|
||||
|
||||
% Init
|
||||
a = alpha;
|
||||
b = beta;
|
||||
n = 0;
|
||||
|
||||
k=1;
|
||||
while b(k) - a(k) > lambda
|
||||
% bisect [a,b]
|
||||
mid = (a(k) + b(k)) / 2;
|
||||
x_1 = mid - epsilon;
|
||||
x_2 = mid + epsilon;
|
||||
|
||||
% set new search interval
|
||||
k = k + 1;
|
||||
if fun(x_1) < fun(x_2)
|
||||
a(k) = a(k-1);
|
||||
b(k) = x_2;
|
||||
else
|
||||
a(k) = x_1;
|
||||
b(k) = b(k-1);
|
||||
end
|
||||
end
|
||||
|
||||
end
|
||||
@@ -0,0 +1,32 @@
|
||||
function [gamma] = gamma_armijo(f, grad_f, dk, xk)
|
||||
% Calculates the best step based on amijo method
|
||||
%
|
||||
% f(xk+ γk*dk) ≤ f(xk) + σ * γk * dk^T * ∇f(xk)
|
||||
% γk = β*γk_0
|
||||
%
|
||||
% f: Objective function
|
||||
% grad_fun: Gradient function of f
|
||||
% dk: Current value of selected direction -∇f or -inv{H}*∇f or -inv{H + lI}*∇f
|
||||
% xk: Current point (x,y)
|
||||
|
||||
|
||||
% beta: beta factor in [0.1, 0.5]
|
||||
% signam: sigma factor in (0, 0.1]
|
||||
global amijo_beta
|
||||
global amijo_sigma
|
||||
|
||||
gf = grad_f(xk);
|
||||
gamma = 1; % Start with a step size of 1
|
||||
|
||||
% Perform Armijo line search
|
||||
while f(xk + gamma * dk) > f(xk) + amijo_sigma * gamma * dk * gf
|
||||
%while f(xk(1) + gamma * dk(1), xk(2) + gamma * dk(2)) > ...
|
||||
% f(xk(1), xk(2)) + amijo_sigma * gamma * norm(dk)^2
|
||||
gamma = amijo_beta * gamma; % Reduce step size
|
||||
if gamma < 1e-12 % Safeguard to prevent infinite reduction
|
||||
warning('Armijo step size became too small.');
|
||||
break;
|
||||
end
|
||||
end
|
||||
|
||||
end
|
||||
@@ -0,0 +1,12 @@
|
||||
function [gamma] = gamma_fixed(~, ~, ~, ~)
|
||||
% Return a fixed step
|
||||
%
|
||||
% This is for completion and code symmetry.
|
||||
%
|
||||
|
||||
global gamma_fixed_step
|
||||
|
||||
% Perform line search
|
||||
gamma = gamma_fixed_step;
|
||||
|
||||
end
|
||||
@@ -0,0 +1,24 @@
|
||||
function [gamma] = gamma_minimized(f, ~, dk, xk)
|
||||
% Calculates the step based on minimizing f(xk− γk*dk)
|
||||
%
|
||||
%
|
||||
% f: Objective function
|
||||
% ~: Gradient function of f - Not used
|
||||
% dk: Current value of selected direction -∇f or -inv{H}*∇f or -inv{H + lI}*∇f
|
||||
% xk: Current point (x,y)
|
||||
|
||||
% Define the line search function fmin(g) = f(xk - g * dk)
|
||||
fmin = @(g) f(xk) + g * dk;
|
||||
|
||||
% find g that minimizes fmin
|
||||
e = 0.0001;
|
||||
l = 0.001;
|
||||
[a,b,k,~] = fmin_bisection(fmin, 0.0001, 1, e, l); % g in (0, 1]
|
||||
gamma = 0.5*(a(k) + b(k));
|
||||
|
||||
% Define the line search function fmin(g) = f(xk - g * dk)
|
||||
%fmin = @(g) f(xk(1) - gamma * dk(1), xk(2) - gamma * dk(2));
|
||||
|
||||
% find g that minimizes fmin
|
||||
%gamma = fminbnd(g, 0, 1);
|
||||
end
|
||||
@@ -0,0 +1,39 @@
|
||||
function [x_vals, f_vals, k] = method_SteepDesc(f, grad_f, xk, tol, max_iter, mode)
|
||||
% f: Objective function
|
||||
% grad_f: Gradient of the function
|
||||
% xk: Initial point [x0; y0]
|
||||
% tol: Tolerance for stopping criterion
|
||||
% max_iter: Maximum number of iterations
|
||||
|
||||
% x_vals: Vector with the (x,y) values until minimum
|
||||
% f_vals: Vector with f(x,y) values until minimum
|
||||
% k: Number of iterations
|
||||
|
||||
if strcmp(mode, 'armijo') == 1
|
||||
gamma_f = @(f, grad_f, dk, xk) gamma_armijo(f, grad_f, dk, xk);
|
||||
elseif strcmp(mode, 'minimized') == 1
|
||||
gamma_f = @(f, grad_f, dk, xk) gamma_minimized(f, grad_f, dk, xk);
|
||||
else % mode == 'fixed'
|
||||
gamma_f = @(f, grad_f, dk, xk) gamma_fixed(f, grad_f, dk, xk);
|
||||
end
|
||||
|
||||
% Storage for iterations, begin with the first point
|
||||
x_vals = xk;
|
||||
f_vals = f(xk);
|
||||
for k = 1:max_iter
|
||||
% Check for convergence
|
||||
if norm(grad_f(xk)) < tol
|
||||
break;
|
||||
end
|
||||
dk = - grad_f(xk); % Steepset descent direction
|
||||
|
||||
gk = gamma_f(f, grad_f, dk, xk); % Calculate gamma
|
||||
|
||||
x_next = xk + gk * dk; % Update step
|
||||
f_next = f(x_next);
|
||||
|
||||
xk = x_next; % Update point
|
||||
x_vals = [x_vals, x_next]; % Store values
|
||||
f_vals = [f_vals, f_next]; % Store function values
|
||||
end
|
||||
end
|
||||
@@ -0,0 +1,50 @@
|
||||
function [x_vals, f_vals, k] = method_SteepDesc_Proj(f, grad_f, xk, sk, limmits, tol, max_iter, mode)
|
||||
% f: Objective function
|
||||
% grad_f: Gradient of the function
|
||||
% xk: Initial point [x0; y0]
|
||||
% sk: Step size (fixed positive scalar)
|
||||
% limits: Bounds of the feasible set for each dimension
|
||||
% tol: Tolerance for stopping criterion
|
||||
% max_iter: Maximum number of iterations
|
||||
|
||||
% x_vals: Vector with the (x,y) values until minimum
|
||||
% f_vals: Vector with f(x,y) values until minimum
|
||||
% k: Number of iterations
|
||||
|
||||
if strcmp(mode, 'armijo') == 1
|
||||
gamma_f = @(f, grad_f, dk, xk) gamma_armijo(f, grad_f, dk, xk);
|
||||
elseif strcmp(mode, 'minimized') == 1
|
||||
gamma_f = @(f, grad_f, dk, xk) gamma_minimized(f, grad_f, dk, xk);
|
||||
else % mode == 'fixed'
|
||||
gamma_f = @(f, grad_f, dk, xk) gamma_fixed(f, grad_f, dk, xk);
|
||||
end
|
||||
|
||||
% Project the first point if needed
|
||||
xk = ProjectionPoint(xk, limmits);
|
||||
|
||||
% Storage for iterations, begin with the first point
|
||||
x_vals = xk;
|
||||
f_vals = f(xk);
|
||||
for k = 1:max_iter
|
||||
% Check for convergence
|
||||
if norm(grad_f(xk)) < tol
|
||||
break;
|
||||
end
|
||||
dk = - grad_f(xk); % Steepset descent direction
|
||||
|
||||
% First calculate xk-bar and project it if nessesary
|
||||
xkbar = xk + sk * dk;
|
||||
xkbar = ProjectionPoint(xkbar, limmits);
|
||||
|
||||
dk = (xkbar - xk); % Steepest descent projection direction
|
||||
|
||||
gk = gamma_f(f, grad_f, dk, xk); % Calculate gamma
|
||||
|
||||
x_next = xk + gk * dk; % Update step
|
||||
f_next = f(x_next);
|
||||
|
||||
xk = x_next; % Update point
|
||||
x_vals = [x_vals, x_next]; % Store values
|
||||
f_vals = [f_vals, f_next]; % Store function values
|
||||
end
|
||||
end
|
||||
@@ -0,0 +1,42 @@
|
||||
function plot3dFun(fun, x_lim, y_lim, size, plot_title, filename)
|
||||
% 3D plots a function
|
||||
% fun: The function to 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
|
||||
%
|
||||
|
||||
global image_width,
|
||||
global image_height;
|
||||
|
||||
% Generate a grid for x and y
|
||||
x_space = linspace(x_lim(1), x_lim(2), size);
|
||||
y_space = linspace(y_lim(1), y_lim(2), size);
|
||||
[X, Y] = meshgrid(x_space, y_space);
|
||||
|
||||
% Evaluate the function on the grid
|
||||
for i = 1:size
|
||||
for j = 1:size
|
||||
% Pass each [x1; x2] as input to fun
|
||||
Z(i, j) = fun([X(i, j); Y(i, j)]);
|
||||
end
|
||||
end
|
||||
|
||||
% 3D plot
|
||||
figure('Name', 'f(x,y)', 'NumberTitle', 'off');
|
||||
set(gcf, 'Position', [100, 100, image_width, image_height]); % Set the figure size
|
||||
surf(X, Y, Z);
|
||||
|
||||
% Customize the plot
|
||||
xlabel('x1'); % Label for x-axis
|
||||
ylabel('x2'); % Label for y-axis
|
||||
zlabel('f(x, y)'); % Label for z-axis
|
||||
title(plot_title, 'Interpreter', 'latex', 'FontSize', 16); % Title of the plot
|
||||
colorbar;
|
||||
|
||||
% save the figure
|
||||
if strcmp(filename, '') == 0
|
||||
print(gcf, filename, '-dpng', '-r300');
|
||||
end
|
||||
end
|
||||
@@ -0,0 +1,41 @@
|
||||
function plotContour(fun, x_lim, y_lim, size, plot_title, filename)
|
||||
% plot the contour of a function
|
||||
% fun: The function to 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
|
||||
%
|
||||
|
||||
global image_width,
|
||||
global image_height;
|
||||
|
||||
% Generate a grid for x and y
|
||||
x_space = linspace(x_lim(1), x_lim(2), size);
|
||||
y_space = linspace(y_lim(1), y_lim(2), size);
|
||||
[X, Y] = meshgrid(x_space, y_space);
|
||||
|
||||
% Evaluate the function on the grid
|
||||
for i = 1:size
|
||||
for j = 1:size
|
||||
% Pass each [x1; x2] as input to fun
|
||||
Z(i, j) = fun([X(i, j); Y(i, j)]);
|
||||
end
|
||||
end
|
||||
|
||||
% Contour
|
||||
figure('Name', 'Contours of f(x1,x2)', 'NumberTitle', 'off');
|
||||
set(gcf, 'Position', [100, 100, image_width, image_height]); % Set the figure size
|
||||
contour(X, Y, Z);
|
||||
|
||||
% Customize the plot
|
||||
xlabel('x1'); % Label for x-axis
|
||||
ylabel('x2'); % Label for y-axis
|
||||
title(plot_title, 'Interpreter', 'latex', 'FontSize', 16); % Title of the plot
|
||||
colorbar;
|
||||
|
||||
% save the figure
|
||||
if strcmp(filename, '') == 0
|
||||
print(gcf, filename, '-dpng', '-r300');
|
||||
end
|
||||
end
|
||||
@@ -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
|
||||
@@ -0,0 +1,28 @@
|
||||
function plotItersOverGamma(gamma, iterations, plot_title, filename)
|
||||
% 3D plots a function
|
||||
% fun: 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;
|
||||
|
||||
figure('Name', 'Iterations_over_gamma', 'NumberTitle', 'off');
|
||||
set(gcf, 'Position', [100, 100, image_width, image_height]); % Set the figure size
|
||||
plot(gamma, iterations, '*r', 'LineWidth', 2);
|
||||
|
||||
% Customize the plot
|
||||
title(plot_title, 'Interpreter', 'latex', 'FontSize', 16); % Title of the plot
|
||||
xlabel('\gamma') ;
|
||||
ylabel('Iterations');
|
||||
|
||||
% save the figure
|
||||
if strcmp(filename, '') == 0
|
||||
print(gcf, filename, '-dpng', '-r300');
|
||||
end
|
||||
end
|
||||
@@ -0,0 +1,48 @@
|
||||
function plotPointsOverContour(points, contour_fun, x_lim, y_lim, size, 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;
|
||||
|
||||
% Generate a grid for x and y
|
||||
x_space = linspace(x_lim(1), x_lim(2), size);
|
||||
y_space = linspace(y_lim(1), y_lim(2), size);
|
||||
[X, Y] = meshgrid(x_space, y_space);
|
||||
|
||||
% Evaluate the function on the grid
|
||||
for i = 1:size
|
||||
for j = 1:size
|
||||
% Pass each [x1; x2] as input to fun
|
||||
Z(i, j) = contour_fun([X(i, j); Y(i, j)]);
|
||||
end
|
||||
end
|
||||
|
||||
% 2D plot
|
||||
figure('Name', '(x1,x2) convergence', 'NumberTitle', 'off');
|
||||
set(gcf, 'Position', [100, 100, image_width, image_height]); % Set the figure size
|
||||
plot(points(1, :), points(2, :), '-or');
|
||||
hold on
|
||||
contour(X, Y, Z);
|
||||
% Customize the plot
|
||||
xlim(x_lim);
|
||||
ylim(y_lim);
|
||||
xlabel('x1'); % Label for x-axis
|
||||
ylabel('x2'); % Label for y-axis
|
||||
grid on
|
||||
|
||||
title(plot_title, 'Interpreter', 'latex', 'FontSize', 16); % Title of the plot
|
||||
colorbar;
|
||||
|
||||
% save the figure
|
||||
if strcmp(filename, '') == 0
|
||||
print(gcf, filename, '-dpng', '-r300');
|
||||
end
|
||||
end
|
||||
Binary file not shown.
Binary file not shown.
Binary file not shown.
@@ -0,0 +1,332 @@
|
||||
%
|
||||
% Optimization Techniques Work 2 report
|
||||
%
|
||||
% authors:
|
||||
% Χρήστος Χουτουρίδης ΑΕΜ 8997
|
||||
% cchoutou@ece.auth.gr
|
||||
|
||||
|
||||
\documentclass[a4paper, 11pt]{AUTHReport}
|
||||
|
||||
% Document configuration
|
||||
\AuthorName{Χρήστος Χουτουρίδης}
|
||||
\AuthorMail{cchoutou@ece.auth.gr}
|
||||
\AuthorAEM{8997}
|
||||
|
||||
% \CoAuthorName{CoAuthor Name}
|
||||
% \CoAuthorMail{CoAuthor Mail}
|
||||
% \CoAuthorAEM{AEM}
|
||||
|
||||
% \WorkGroup{Ομάδα Χ}
|
||||
|
||||
\DocTitle{2η Εργαστηριακή Άσκηση}
|
||||
\DocSubTitle{Ελαχιστοποίηση συναρτήσεων πολλών μεταβλητών χωρίς περιορισμούς με χρήση παραγώγων}
|
||||
|
||||
\Department{Τμήμα ΗΜΜΥ. Τομέας Ηλεκτρονικής}
|
||||
\ClassName{Τεχνικές Βελτιστοποίησης}
|
||||
|
||||
\InstructorName{Γ. Ροβιθάκης}
|
||||
\InstructorMail{rovithak@auth.gr}
|
||||
\CurrentDate{\today}
|
||||
|
||||
\usepackage{capt-of}
|
||||
\usepackage{enumitem}
|
||||
\usepackage{tabularx}
|
||||
\usepackage{array}
|
||||
\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}
|
||||
|
||||
\InsertTitle
|
||||
|
||||
%\tableofcontents
|
||||
|
||||
\sloppy
|
||||
|
||||
\section{Εισαγωγή}
|
||||
Η παρούσα εργασία αφορά το πρόβλημα της ελαχιστοποίησης μιας δοσμένης συνάρτησης πολλών μεταβλητών $f: \mathbb{R}^n \rightarrow \mathbb{R}$ χωρίς περιορισμούς.
|
||||
Για το σκοπό αυτό κάνουμε χρήση τριών μεθόδων.
|
||||
Της μεθόδου μέγιστης καθόδου (Steepest Descent), της μεθόδου Newton, και της Levenberg-Marquardt.
|
||||
Ακόμα για κάθε μία από αυτές θα υλοποιήσουμε τρεις διαφορετικές τεχνικές υπολογισμού βήματος.
|
||||
|
||||
\section{Παραδοτέα}
|
||||
Τα παραδοτέα της εργασίας αποτελούνται από:
|
||||
\begin{itemize}
|
||||
\item Την παρούσα αναφορά.
|
||||
\item Τον κατάλογο \textbf{scripts/}, που περιέχει τον κώδικα της MATLAB.
|
||||
\item Το \href{https://git.hoo2.net/hoo2/OptimizationTechniques/src/branch/master/Work2}{σύνδεσμο} με το αποθετήριο που περιέχει όλο το project με τον κώδικα της MATLAB, της αναφοράς και τα παραδοτέα.
|
||||
\end{itemize}
|
||||
|
||||
\section{Προγραμματιστική προσέγγιση}
|
||||
Για τον προγραμματισμό και εκτέλεση των μεθόδων της παρούσας εργασίας έγινε χρήση της MATLAB.
|
||||
Στον κατάλογο \textbf{scripts}, περιέχονται όλες οι μέθοδοι και οι τεχνικές υπολογισμού βημάτων με τη μορφή συναρτήσεων καθώς και scripts που τις καλούν.
|
||||
Για κάθε μία μέθοδο (ένα θέμα της εργασίας), υπάρχει το αντίστοιχο script που περιέχει τους υπολογισμούς, τις κλήσεις των μεθόδων και τη δημιουργία των διαγραμμάτων.
|
||||
Για το πρώτο θέμα το αρχείο Script\_1\_Plots.m για το δεύτερο το Script\_2\_Steepest\_descent.m και ούτω καθεξής.
|
||||
Στην παρούσα εργασία η υλοποίηση του κώδικα ακολουθεί την τεχνική της προηγούμενης εργασίας και “ομαδοποιεί” αρκετές λειτουργίες.
|
||||
Πιο συγκεκριμένα.
|
||||
|
||||
\subsection{Κλήση μεθόδων επιλογής βήματος $\gamma_k$}
|
||||
\label{subsec:polymorphic-calls}
|
||||
Δεδομένου ότι οι μέθοδοι θα πρέπει να καλεστούν και εκτελεστούν με παραπάνω από μία τεχνική επιλογής βήματος $\gamma_k$, δημιουργήσαμε εσωτερικά της κάθε μεθόδου ένα κοινό interface για τις μεθόδους επιλογής βήματος.
|
||||
Αυτό έχει τη μορφή: \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 και με βάση αυτή γίνεται η τελική επιλογή.
|
||||
Έτσι για παράδειγμα η κλήση \textit{method\_newtown(f, gradf, [0, 0], 0.001, 1000, 'armijo')} υλοποιεί τη μέθοδο newton χρησιμοποιώντας τη μέθοδο Armijo για επιλογή βήματος ενώ η \textit{method\_newtown(f, gradf, [0, 0], 0.001, 1000, 'minimized')}, χρησιμοποιεί βήμα που ελαχιστοποιεί την $f(x_k + \gamma_k d_k)$.
|
||||
|
||||
\subsection{Symbolic expression functions}
|
||||
Μία ακόμη προγραμματιστική τεχνική που ακολουθήθηκε είναι η χρήση \textbf{symbolic expression} για την αναπαράσταση των διαφορετικών αντικειμενικών συναρτήσεων.
|
||||
Ο λόγος που επιλέχθηκε είναι η \textbf{δυνατότητα εξαγωγής ενός symbolic expression που αναπαριστά την κλίση $\nabla f$ και τον Εσσιανό $\nabla^2f$ μιας συνάρτησης} από την MATLAB, κάνοντας χρήση των εντολών \textit{gradient()} και \textit{hessian()}.
|
||||
Αν αντίθετα χρησιμοποιούσαμε απλές συναρτήσεις, πολυώνυμα ή lambdas για την αναπαράσταση των αντικειμενικών συναρτήσεων, τότε για τον υπολογισμό της κλίσης και του Εσσιανού θα έπρεπε:
|
||||
\begin{itemize}
|
||||
\item Είτε να υπολογίζαμε αριθμητικά τις παραγώγους gradient και hessian μέσα στις μεθόδους, κάτι που θα εισήγαγε \textit{\textbf{αχρείαστο αριθμητικό σφάλμα}}.
|
||||
\item Είτε να κάναμε χρήση δύο επιπλέων συναρτήσεων (ή πολυωνύμων) για την αναπαράσταση τους, κάτι που ουσιαστικά θα δημιουργούσε \textit{\textbf{πλεονασμό πληροφορίας εισόδου}} και άρα μεγαλύτερη πιθανότητα να κάνουμε λάθος.
|
||||
\end{itemize}
|
||||
Η αναπαράσταση όμως με χρήση symbolic expression είναι πιο “βαριά” όταν χρειάζεται να υπολογίσουμε την τιμή μιας συνάρτησης σε κάποιο σημείο (subs(expr, number)).
|
||||
Αυτό είναι κάτι που χρειάζεται εκτενώς στον κώδικά μας.
|
||||
Για το λόγο αυτό, ενώ η συνάρτηση δίνεται ως symbolic expression, μέσω αυτής υπολογίζονται αυτόματα η κλίση, ο Εσσιανός αλλά και οι “κανονικές” συναρτήσεις MATLAB που τις υλοποιούν.
|
||||
Έτσι έχουμε την ακριβή αναπαράσταση της κλίσης και του Εσσιανού ως συναρτήσεις χωρίς να πληρώνουμε το κόστος της subs().
|
||||
|
||||
\section{Απεικόνιση της συνάρτησης - Θέμα 1}
|
||||
Η συνάρτηση με την οποία ασχολούμαστε στην παρούσα εργασία είναι η:
|
||||
\boldmath
|
||||
\begin{equation}
|
||||
f(x,y) = x^5 \cdot e^{-x^2 - y^2}
|
||||
\end{equation}
|
||||
\label{eq:ObjectiveFunction}
|
||||
\unboldmath
|
||||
Στο παρακάτω σχήμα \ref{fig:plot3dFunction} φαίνεται η τρισδιάστατη απεικόνιση της συνάρτησης.
|
||||
\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/Plot_Contour.png}{Ισοβαρείς της f}
|
||||
|
||||
Από το παραπάνω σχήμα \ref{fig:plotContour} φαίνεται ότι το ελάχιστο της f βρίσκεται στο αρνητικό ημιεπίπεδο των χ, κοντά στο $y = 0$
|
||||
\par
|
||||
Τα διαγράμματα για τη μέθοδο δημιουργούνται εκτελώντας το αρχείο \textbf{Script\_1\_Plots.m}
|
||||
|
||||
\section{Τεχνικές υπολογισμού βήματος}
|
||||
Πριν προχωρήσουμε στα επόμενα θέματα της εργασίας και στην ανάλυση των μεθόδων υπολογισμού του ελάχιστου, θέλουμε να αναφερθούμε στις διαφορετικές τεχνικές επιλογής βήματος $\gamma_k$ και ειδικότερα για αυτή της ελαχιστοποίησης της $f(x_k + \gamma_k d_k)$ και την Armijo.
|
||||
|
||||
\subsection{Ελαχιστοποίηση της $f(x_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)\} \]
|
||||
Η κατεύθυνση $d_k$ μπορεί να είναι:
|
||||
\begin{itemize}
|
||||
\item Η αρνητική κλίση $-\nabla f(x_k)$ (Steepest Descent).
|
||||
\item Η Newton direction $-{H_k}^{-1} \nabla f(x_k)$ (Newton και Levenberg-Marquardt).
|
||||
\end{itemize}
|
||||
\unboldmath
|
||||
Πλεονεκτήματα της μεθόδου είναι η \textbf{ακρίβεια}, καθώς το βήμα $\gamma_k$ υπολογίζεται με βέλτιστο τρόπο για τη συγκεκριμένη κατεύθυνση και η \textbf{γρήγορη σύγκλιση}, ειδικά σε προβλήματα όπου η $f(x)$ είναι καμπυλωτή με $\gamma_k$ επακριβώς ορισμένο.
|
||||
Στα μειονεκτήματα μπορούμε να αναφέρουμε το \textbf{υπολογιστικό κόστος}, καθώς ο υπολογισμός του $\gamma_k$ απαιτεί πολλαπλές αξιολογήσεις της $f(x)$.
|
||||
|
||||
\subsection{Armijo rule}
|
||||
Η 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) \]
|
||||
Όπου $\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{Μέθοδος Μέγιστης Καθόδου - Θέμα 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$ να είναι υπολογίσιμη.
|
||||
Επίσης για την εφαρμογή της μεθόδου, για το αρχικό σημείο θα πρέπει \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.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:StDes_Iter_o_gamma_2}{../scripts/figures/StDes_Iter_o_gamma_2.png}{Αριθμός επαναλήψεων για διαφορετικές τιμές $\gamma_k$ [Μέγιστη Κάθοδος].}
|
||||
|
||||
Στο παραπάνω σχήμα \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.1353$ και το διάνυσμα της κλίσης: $\nabla f(0,0) = \begin{bmatrix} 0.4060 \\ 0.2707 \end{bmatrix}$, επομένως μπορούμε να εφαρμόσουμε τη μέθοδο.
|
||||
|
||||
\par
|
||||
\underline{Σταθερό βήμα} \\
|
||||
Για σταθερό βήμα εκτελέσαμε διαδοχικά τη μέθοδο \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].}
|
||||
Για υπολογισμό βήματος ελαχιστοποιώντας τη $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 δεν συγκλίνει, καθώς και αυτή εγκλωβίζεται στο δεξιό ημιεπίπεδο.
|
||||
|
||||
\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(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.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.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 (αρχείο: \textbf{method\_lev\_mar.m}).
|
||||
Πρόκειται για μια τροποποιημένη έκδοση της μεθόδου Newton, η οποία εισάγει έναν παράγοντα απόσβεσης για τη σταθεροποίηση όταν ο εσσιανός δεν είναι θετικά ορισμένος.
|
||||
\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}
|
||||
|
||||
\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 = (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$].}
|
||||
|
||||
\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}
|
||||
@@ -0,0 +1,33 @@
|
||||
% Given environment
|
||||
|
||||
clear;
|
||||
% Setup the function under test
|
||||
syms x y;
|
||||
fexpr = x^5 * exp(-x^2 - y^2);
|
||||
title_fun = "$f(x,y) = x^5 \cdot e^{-x^2 - y^2}$";
|
||||
|
||||
% Calculate the gradient and Hessian
|
||||
grad_fexpr = gradient(fexpr, [x, y]); % Gradient of f
|
||||
hessian_fexpr = hessian(fexpr, [x, y]); % Hessian of f
|
||||
|
||||
% Convert symbolic expressions to MATLAB functions
|
||||
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])
|
||||
|
||||
%fixed step size globals
|
||||
global gamma_fixed_step
|
||||
|
||||
global image_width,
|
||||
global image_height;
|
||||
|
||||
image_width = 960;
|
||||
image_height = 640;
|
||||
@@ -0,0 +1,14 @@
|
||||
% Define environment (functions, gradients etc...)
|
||||
GivenEnv
|
||||
|
||||
%
|
||||
% We plot the function in the domain of x,y in [-3, 3].
|
||||
% We also plot the contour in order to get a sense of the min and maximum
|
||||
% points in the x-y plane
|
||||
%
|
||||
|
||||
% 3d plot the function
|
||||
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/Plot_Contour.png");
|
||||
@@ -0,0 +1,116 @@
|
||||
% Define environment (functions, gradients etc...)
|
||||
GivenEnv
|
||||
|
||||
% Define parameters
|
||||
max_iter = 300; % Maximum iterations
|
||||
tol = 1e-4; % Tolerance
|
||||
|
||||
% Point x0 = (0, 0)
|
||||
% =========================================================================
|
||||
point = 1;
|
||||
x0 = [0, 0];
|
||||
f = fun(x0(1), x0(2));
|
||||
gf = grad_fun(x0(1), x0(2));
|
||||
hf = hessian_fun(x0(1), x0(2));
|
||||
fprintf('Initial point (%d, %d), f = %f, grad = [%f;%f], hessian = [%f %f ; %f %f]. Can NOT use method\n', x0, f, gf, hf);
|
||||
disp(' ');
|
||||
|
||||
|
||||
% Point x0 = (-1, 1)
|
||||
% =========================================================================
|
||||
point = 2;
|
||||
x0 = [-1, 1];
|
||||
point_str = "[" + x0(1) + ", " + x0(2) + "]";
|
||||
|
||||
f = fun(-1, 1);
|
||||
gf = grad_fun(x0(1), x0(2));
|
||||
hf = hessian_fun(x0(1), x0(2));
|
||||
fprintf('Initial point (%d, %d), f = %f, grad = [%f;%f], hessian = [%f %f ; %f %f]. Can use method\n', x0, f, gf, hf);
|
||||
|
||||
% Find the best fixed gamma
|
||||
k = zeros(100, 1);
|
||||
j = 1;
|
||||
n = linspace(0.1, 1.5, 100);
|
||||
for g = n
|
||||
gamma_fixed_step = g;
|
||||
[~, ~, k(j)] = method_steepest_descent(fun, grad_fun, x0, tol, max_iter, 'fixed');
|
||||
j = j + 1;
|
||||
end
|
||||
plotItersOverGamma(n, k, "Iteration for different $\gamma$ values", "figures/StDes_Iter_o_gamma_" + point + ".png");
|
||||
|
||||
[~, j] = min(k);
|
||||
gamma_fixed_step = n(j);
|
||||
|
||||
[x_fixed, f_fixed, kk] = method_steepest_descent(fun, grad_fun, x0, tol, max_iter, 'fixed');
|
||||
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, [-2, 0], [-2, 2], 100, point_str + ": Steepest descent $\gamma$ = " + gamma_fixed_step, "figures/StDes_fixed_" + point + ".png");
|
||||
|
||||
|
||||
% Minimized f
|
||||
[x_minimized, f_minimized, kk] = method_steepest_descent(fun, grad_fun, 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, [-2, 0], [-2, 2], 100, point_str + ": Steepest descent minimized $f(x_k + \gamma_kd_k)$", "figures/StDes_minimized_" + point + ".png");
|
||||
|
||||
|
||||
% Armijo Rule
|
||||
|
||||
% Methods tuning
|
||||
amijo_beta = 0.4; % typical range: [0.1, 0.8]
|
||||
amijo_sigma = 0.1; % typical range: [0.01, 0.3]
|
||||
|
||||
[x_armijo, f_armijo, kk] = method_steepest_descent(fun, grad_fun, 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, [-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;
|
||||
x0 = [1, -1];
|
||||
point_str = "[" + x0(1) + ", " + x0(2) + "]";
|
||||
|
||||
f = fun(-1, 1);
|
||||
gf = grad_fun(x0(1), x0(2));
|
||||
hf = hessian_fun(x0(1), x0(2));
|
||||
fprintf('Initial point (%d, %d), f = %f, grad = [%f;%f], hessian = [%f %f ; %f %f]. Can use method\n', x0, f, gf, hf);
|
||||
|
||||
% Find the best fixed gamma
|
||||
k = zeros(100, 1);
|
||||
j = 1;
|
||||
n = linspace(0.1, 1, 100);
|
||||
for g = n
|
||||
gamma_fixed_step = g;
|
||||
[~, ~, k(j)] = method_steepest_descent(fun, grad_fun, x0, tol, max_iter, 'fixed');
|
||||
j = j + 1;
|
||||
end
|
||||
%if min(k) == max_iter
|
||||
% fprintf('Fixed step: Initial point (%d, %d). Can NOT use method\n', x0);
|
||||
%end
|
||||
%plotItersOverGamma(n, k, "Iteration for different $\gamma$ values", "figures/StDes_Iter_o_gamma_" + point + ".png");
|
||||
|
||||
%gamma_fixed_step = 1;
|
||||
[x_fixed, f_fixed, kk] = method_steepest_descent(fun, grad_fun, x0, tol, max_iter, 'fixed');
|
||||
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, [-1, 2], [-2, 2], 100, point_str + ": Steepest descent $\gamma$ = " + gamma_fixed_step, "figures/StDes_fixed_" + point + ".png");
|
||||
|
||||
|
||||
% Minimized f
|
||||
[x_minimized, f_minimized, kk] = method_steepest_descent(fun, grad_fun, 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, [-2, 2], [-3, 2], 100, point_str + ": Steepest descent minimized $f(x_k + \gamma_kd_k)$", "figures/StDes_minimized_" + point + ".png");
|
||||
|
||||
|
||||
% Armijo Rule
|
||||
|
||||
% Methods tuning
|
||||
amijo_beta = 0.4; % typical range: [0.1, 0.8]
|
||||
amijo_sigma = 0.1; % typical range: [0.01, 0.3]
|
||||
|
||||
[x_armijo, f_armijo, kk] = method_steepest_descent(fun, grad_fun, 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, [-1, 2], [-2, 2], 100, point_str + ": Steepest descent Armijo method", "figures/StDes_armijo_" + point + ".png");
|
||||
|
||||
|
||||
@@ -0,0 +1,42 @@
|
||||
% Define environment (functions, gradients etc...)
|
||||
GivenEnv
|
||||
|
||||
|
||||
% Point x0 = (0, 0)
|
||||
% =========================================================================
|
||||
point = 1;
|
||||
x0 = [0, 0];
|
||||
f = fun(x0(1), x0(2));
|
||||
gf = grad_fun(x0(1), x0(2));
|
||||
hf = hessian_fun(x0(1), x0(2));
|
||||
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);
|
||||
disp(' ');
|
||||
|
||||
|
||||
% Point x0 = (-1, 1)
|
||||
% =========================================================================
|
||||
point = 2;
|
||||
x0 = [-1, 1];
|
||||
point_str = "[" + x0(1) + ", " + x0(2) + "]";
|
||||
|
||||
f = fun(-1, 1);
|
||||
gf = grad_fun(x0(1), x0(2));
|
||||
hf = hessian_fun(x0(1), x0(2));
|
||||
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);
|
||||
disp(' ');
|
||||
|
||||
% Point x0 = (1, -1)
|
||||
% =========================================================================
|
||||
point = 3;
|
||||
x0 = [1, -1];
|
||||
point_str = "[" + x0(1) + ", " + x0(2) + "]";
|
||||
|
||||
f = fun(-1, 1);
|
||||
gf = grad_fun(x0(1), x0(2));
|
||||
hf = hessian_fun(x0(1), x0(2));
|
||||
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);
|
||||
|
||||
|
||||
@@ -0,0 +1,119 @@
|
||||
% Define environment (functions, gradients etc...)
|
||||
GivenEnv
|
||||
|
||||
% Define parameters
|
||||
max_iter = 300; % Maximum iterations
|
||||
tol = 1e-4; % Tolerance
|
||||
|
||||
|
||||
% Point x0 = (0, 0)
|
||||
% =========================================================================
|
||||
point = 1;
|
||||
x0 = [0, 0];
|
||||
f = fun(x0(1), x0(2));
|
||||
gf = grad_fun(x0(1), x0(2));
|
||||
hf = hessian_fun(x0(1), x0(2));
|
||||
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);
|
||||
disp(' ');
|
||||
|
||||
|
||||
% Point x0 = (-1, 1)
|
||||
% =========================================================================
|
||||
point = 2;
|
||||
x0 = [-1, 1];
|
||||
point_str = "[" + x0(1) + ", " + x0(2) + "]";
|
||||
|
||||
f = fun(-1, 1);
|
||||
gf = grad_fun(x0(1), x0(2));
|
||||
hf = hessian_fun(x0(1), x0(2));
|
||||
ev = eig(hf);
|
||||
fprintf('Initial point (%d, %d), f = %f, grad = [%f;%f], hessian = [%f %f ; %f %f]. Eigenvalues= [%f, %f], Can use method\n', x0, f, gf, hf, ev);
|
||||
|
||||
|
||||
% Find the best fixed gamma
|
||||
k = zeros(100, 1);
|
||||
j = 1;
|
||||
n = linspace(0.1, 1.5, 100);
|
||||
for g = n
|
||||
gamma_fixed_step = g;
|
||||
[x, f, k(j)] = method_lev_mar(fun, grad_fun, hessian_fun, 0.3, x0, tol, max_iter, 'fixed');
|
||||
if ~(x(end, 1) < -1.57 && x(end, 1) > -1.59 && x(end, 2) < 0.01 && x(end,2) > -0.01 && f(end) < -0.8 && f(end) > -0.82)
|
||||
k(j) = 300;
|
||||
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);
|
||||
|
||||
[x_fixed, f_fixed, kk] = method_lev_mar(fun, grad_fun, hessian_fun, 0.3, x0, tol, max_iter, 'fixed');
|
||||
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_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
|
||||
|
||||
% Methods tuning
|
||||
amijo_beta = 0.4; % typical range: [0.1, 0.8]
|
||||
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/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)
|
||||
% =========================================================================
|
||||
point = 3;
|
||||
x0 = [1, -1];
|
||||
point_str = "[" + x0(1) + ", " + x0(2) + "]";
|
||||
|
||||
f = fun(-1, 1);
|
||||
gf = grad_fun(x0(1), x0(2));
|
||||
hf = hessian_fun(x0(1), x0(2));
|
||||
ev = eig(hf);
|
||||
fprintf('Initial point (%d, %d), f = %f, grad = [%f;%f], hessian = [%f %f ; %f %f]. Eigenvalues= [%f, %f], Can use method\n', x0, f, gf, hf, ev);
|
||||
|
||||
|
||||
% Find the best fixed gamma
|
||||
k = zeros(100, 1);
|
||||
j = 1;
|
||||
n = linspace(0.1, 1.5, 100);
|
||||
for g = n
|
||||
gamma_fixed_step = g;
|
||||
[x, f, k(j)] = method_lev_mar(fun, grad_fun, hessian_fun, 0.3, x0, tol, max_iter, 'fixed');
|
||||
if ~(x(end, 1) < -1.57 && x(end, 1) > -1.59 && x(end, 2) < 0.01 && x(end,2) > -0.01 && f(end) < -0.8 && f(end) > -0.82)
|
||||
k(j) = 300;
|
||||
end
|
||||
j = j + 1;
|
||||
end
|
||||
|
||||
[~, j] = min(k);
|
||||
gamma_fixed_step = n(j);
|
||||
|
||||
[x_fixed, f_fixed, kk] = method_lev_mar(fun, grad_fun, hessian_fun, 0.3, x0, tol, max_iter, 'fixed');
|
||||
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, 2], [-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, 2], [-2, 2], 100, point_str + ": Levenberg-Marquardt minimized $f(x_k + \gamma_kd_k)$", "figures/LevMar_minimized_" + point + ".png");
|
||||
|
||||
% Armijo Rule
|
||||
|
||||
% Methods tuning
|
||||
amijo_beta = 0.4; % typical range: [0.1, 0.8]
|
||||
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/LevMar_armijo_" + point + ".png");
|
||||
disp(' ');
|
||||
Binary file not shown.
|
After Width: | Height: | Size: 57 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 98 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 86 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 72 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 100 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 86 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 103 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 93 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 73 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 256 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 55 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 93 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 86 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 80 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 94 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 84 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 98 KiB |
Binary file not shown.
|
After Width: | Height: | Size: 92 KiB |
@@ -0,0 +1,49 @@
|
||||
function [a, b, k, n] = fmin_bisection(fun, alpha, beta, epsilon, lambda)
|
||||
% Bisection method for finding the local minimum of a function.
|
||||
%
|
||||
% fun: The objective function
|
||||
% alpha: (number) The starting point of the interval in which we seek
|
||||
% for minimum
|
||||
% beta: (number) The ending point of the interval in which we seek
|
||||
% for minimum
|
||||
% epsilon: (number) The epsilon value (distance from midpoint)
|
||||
% lambda: (number) The lambda value (accuracy)
|
||||
%
|
||||
% return:
|
||||
% a: (vector) Starting points of the interval for each iteration
|
||||
% b: (vector) Ending points of the interval for each iteration
|
||||
% k: (number) The number of iterations
|
||||
% n: (number) The calls of objective function fun_expr
|
||||
%
|
||||
|
||||
|
||||
|
||||
% Error checking
|
||||
if alpha > beta || 2*epsilon >= lambda || lambda <= 0
|
||||
error ('Input criteria not met')
|
||||
end
|
||||
|
||||
% Init
|
||||
a = alpha;
|
||||
b = beta;
|
||||
n = 0;
|
||||
|
||||
k=1;
|
||||
while b(k) - a(k) > lambda
|
||||
% bisect [a,b]
|
||||
mid = (a(k) + b(k)) / 2;
|
||||
x_1 = mid - epsilon;
|
||||
x_2 = mid + epsilon;
|
||||
|
||||
% set new search interval
|
||||
k = k + 1;
|
||||
if fun(x_1) < fun(x_2)
|
||||
a(k) = a(k-1);
|
||||
b(k) = x_2;
|
||||
else
|
||||
a(k) = x_1;
|
||||
b(k) = b(k-1);
|
||||
end
|
||||
end
|
||||
|
||||
end
|
||||
@@ -0,0 +1,33 @@
|
||||
function [gamma] = gamma_armijo(f, grad_f, dk, xk)
|
||||
% Calculates the best step based on amijo method
|
||||
%
|
||||
% f(xk+ γk*dk) ≤ f(xk) + σ * γk * dk^T * ∇f(xk)
|
||||
% γk = β*γk_0
|
||||
%
|
||||
% f: Objective function
|
||||
% grad_fun: Gradient function of f
|
||||
% dk: Current value of selected direction -∇f or -inv{H}*∇f or -inv{H + lI}*∇f
|
||||
% xk: Current point (x,y)
|
||||
|
||||
|
||||
% beta: beta factor in [0.1, 0.5]
|
||||
% signam: sigma factor in (0, 0.1]
|
||||
global amijo_beta
|
||||
global amijo_sigma
|
||||
|
||||
gf = grad_f(xk(1), xk(2));
|
||||
gamma = 1; % Start with a step size of 1
|
||||
|
||||
% Perform Armijo line search
|
||||
while f(xk(1) + gamma * dk(1), xk(2) + gamma * dk(2)) > ...
|
||||
f(xk(1), xk(2)) + amijo_sigma * gamma * dk' * gf
|
||||
%while f(xk(1) + gamma * dk(1), xk(2) + gamma * dk(2)) > ...
|
||||
% f(xk(1), xk(2)) + amijo_sigma * gamma * norm(dk)^2
|
||||
gamma = amijo_beta * gamma; % Reduce step size
|
||||
if gamma < 1e-12 % Safeguard to prevent infinite reduction
|
||||
warning('Armijo step size became too small.');
|
||||
break;
|
||||
end
|
||||
end
|
||||
|
||||
end
|
||||
@@ -0,0 +1,12 @@
|
||||
function [gamma] = gamma_fixed(~, ~, ~, ~)
|
||||
% Return a fixed step
|
||||
%
|
||||
% This is for completion and code symmetry.
|
||||
%
|
||||
|
||||
global gamma_fixed_step
|
||||
|
||||
% Perform line search
|
||||
gamma = gamma_fixed_step;
|
||||
|
||||
end
|
||||
@@ -0,0 +1,24 @@
|
||||
function [gamma] = gamma_minimized(f, ~, dk, xk)
|
||||
% Calculates the step based on minimizing f(xk− γk*dk)
|
||||
%
|
||||
%
|
||||
% f: Objective function
|
||||
% ~: Gradient function of f - Not used
|
||||
% dk: Current value of selected direction -∇f or -inv{H}*∇f or -inv{H + lI}*∇f
|
||||
% xk: Current point (x,y)
|
||||
|
||||
% Define the line search function fmin(g) = f(xk - g * dk)
|
||||
fmin = @(g) f(xk(1) + g*dk(1), xk(2) + g*dk(2));
|
||||
|
||||
% find g that minimizes fmin
|
||||
e = 0.0001;
|
||||
l = 0.001;
|
||||
[a,b,k,~] = fmin_bisection(fmin, 0, 5, e, l);
|
||||
gamma = 0.5*(a(k) + b(k));
|
||||
|
||||
% Define the line search function fmin(g) = f(xk - g * dk)
|
||||
%fmin = @(g) f(xk(1) - gamma * dk(1), xk(2) - gamma * dk(2));
|
||||
|
||||
% find g that minimizes fmin
|
||||
%gamma = fminbnd(g, 0, 1);
|
||||
end
|
||||
@@ -0,0 +1,61 @@
|
||||
function [x_vals, f_vals, k] = method_lev_mar(f, grad_f, hessian_f, e, xk, tol, max_iter, mode)
|
||||
% f: Objective function
|
||||
% grad_f: Gradient of the function
|
||||
% hessian_f: Hessian of the function
|
||||
% 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
|
||||
|
||||
% x_vals: Vector with the (x,y) values until minimum
|
||||
% f_vals: Vector with f(x,y) values until minimum
|
||||
% k: Number of iterations
|
||||
|
||||
|
||||
if strcmp(mode, 'armijo') == 1
|
||||
gamma_f = @(f, grad_f, dk, xk) gamma_armijo(f, grad_f, dk, xk);
|
||||
elseif strcmp(mode, 'minimized') == 1
|
||||
gamma_f = @(f, grad_f, dk, xk) gamma_minimized(f, grad_f, dk, xk);
|
||||
else % mode == 'fixed'
|
||||
gamma_f = @(f, grad_f, dk, xk) gamma_fixed(f, grad_f, dk, xk);
|
||||
end
|
||||
|
||||
x_vals = xk; % Store iterations
|
||||
f_vals = f(xk(1), xk(2));
|
||||
|
||||
for k = 1:max_iter
|
||||
grad = grad_f(xk(1), xk(2));
|
||||
|
||||
% Check for convergence
|
||||
if norm(grad) < tol
|
||||
break;
|
||||
end
|
||||
hess = hessian_f(xk(1), xk(2));
|
||||
|
||||
% 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));
|
||||
if min(eig(hess + mI)) <= 0 % Fail-check
|
||||
warning('Can not normalize hessian matrix.');
|
||||
end
|
||||
end
|
||||
|
||||
% Solve for search direction using Newton's step
|
||||
dk = - inv(hess + mI) * grad;
|
||||
|
||||
% Calculate gamma
|
||||
gk = gamma_f(f, grad_f, dk, xk);
|
||||
|
||||
x_next = xk + gk * dk'; % Update step
|
||||
f_next = f(x_next(1), x_next(2));
|
||||
|
||||
xk = x_next; % Update point
|
||||
x_vals = [x_vals; x_next]; % Store values
|
||||
f_vals = [f_vals; f_next]; % Store function values
|
||||
end
|
||||
end
|
||||
@@ -0,0 +1,47 @@
|
||||
function [x_vals, f_vals, k] = method_newton(f, grad_f, hessian_f, xk, tol, max_iter, mode)
|
||||
% f: Objective function
|
||||
% grad_f: Gradient of the function
|
||||
% hessian_f: Hessian of the function
|
||||
% x0: Initial point [x0, y0]
|
||||
% tol: Tolerance for stopping criterion
|
||||
% max_iter: Maximum number of iterations
|
||||
|
||||
% x_vals: Vector with the (x,y) values until minimum
|
||||
% f_vals: Vector with f(x,y) values until minimum
|
||||
% k: Number of iterations
|
||||
|
||||
|
||||
if strcmp(mode, 'armijo') == 1
|
||||
gamma_f = @(f, grad_f, dk, xk) gamma_armijo(f, grad_f, dk, xk);
|
||||
elseif strcmp(mode, 'minimized') == 1
|
||||
gamma_f = @(f, grad_f, dk, xk) gamma_minimized(f, grad_f, dk, xk);
|
||||
else % mode == 'fixed'
|
||||
gamma_f = @(f, grad_f, dk, xk) gamma_fixed(f, grad_f, dk, xk);
|
||||
end
|
||||
|
||||
x_vals = xk; % Store iterations
|
||||
f_vals = f(xk(1), xk(2));
|
||||
|
||||
for k = 1:max_iter
|
||||
grad = grad_f(xk(1), xk(2));
|
||||
|
||||
% Check for convergence
|
||||
if norm(grad) < tol
|
||||
break;
|
||||
end
|
||||
hess = hessian_f(xk(1), xk(2));
|
||||
|
||||
% Solve for search direction using Newton's step
|
||||
dk = - inv(hess) * grad;
|
||||
|
||||
% Calculate gamma
|
||||
gk = gamma_f(f, grad_f, dk, xk);
|
||||
|
||||
x_next = xk + gk * dk'; % Update step
|
||||
f_next = f(x_next(1), x_next(2));
|
||||
|
||||
xk = x_next; % Update point
|
||||
x_vals = [x_vals; x_next]; % Store values
|
||||
f_vals = [f_vals; f_next]; % Store function values
|
||||
end
|
||||
end
|
||||
@@ -0,0 +1,43 @@
|
||||
function [x_vals, f_vals, k] = method_steepest_descent(f, grad_f, xk, tol, max_iter, mode)
|
||||
% f: Objective function
|
||||
% grad_f: Gradient of the function
|
||||
% xk: Initial point [x0, y0]
|
||||
% tol: Tolerance for stopping criterion
|
||||
% max_iter: Maximum number of iterations
|
||||
|
||||
% x_vals: Vector with the (x,y) values until minimum
|
||||
% f_vals: Vector with f(x,y) values until minimum
|
||||
% k: Number of iterations
|
||||
|
||||
if strcmp(mode, 'armijo') == 1
|
||||
gamma_f = @(f, grad_f, dk, xk) gamma_armijo(f, grad_f, dk, xk);
|
||||
elseif strcmp(mode, 'minimized') == 1
|
||||
gamma_f = @(f, grad_f, dk, xk) gamma_minimized(f, grad_f, dk, xk);
|
||||
else % mode == 'fixed'
|
||||
gamma_f = @(f, grad_f, dk, xk) gamma_fixed(f, grad_f, dk, xk);
|
||||
end
|
||||
|
||||
% Storage for iterations, begin with the first point
|
||||
x_vals = xk;
|
||||
f_vals = f(xk(1), xk(2));
|
||||
|
||||
for k = 1:max_iter
|
||||
grad = grad_f(xk(1), xk(2));
|
||||
|
||||
% Check for convergence
|
||||
if norm(grad) < tol
|
||||
break;
|
||||
end
|
||||
dk = - grad;
|
||||
|
||||
% Calculate gamma
|
||||
gk = gamma_f(f, grad_f, dk, xk);
|
||||
|
||||
x_next = xk + gk * dk'; % Update step
|
||||
f_next = f(x_next(1), x_next(2));
|
||||
|
||||
xk = x_next; % Update point
|
||||
x_vals = [x_vals; x_next]; % Store values
|
||||
f_vals = [f_vals; f_next]; % Store function values
|
||||
end
|
||||
end
|
||||
@@ -0,0 +1,37 @@
|
||||
function plot3dFun(fun, x_lim, y_lim, size, plot_title, filename)
|
||||
% 3D plots a function
|
||||
% fun: The function to 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
|
||||
%
|
||||
|
||||
global image_width,
|
||||
global image_height;
|
||||
|
||||
% Generate a grid for x and y
|
||||
x_space = linspace(x_lim(1), x_lim(2), size);
|
||||
y_space = linspace(y_lim(1), y_lim(2), size);
|
||||
[X, Y] = meshgrid(x_space, y_space);
|
||||
|
||||
% Evaluate the function on the grid
|
||||
Z = fun(X, Y);
|
||||
|
||||
% 3D plot
|
||||
figure('Name', 'f(x,y)', 'NumberTitle', 'off');
|
||||
set(gcf, 'Position', [100, 100, image_width, image_height]); % Set the figure size
|
||||
surf(X, Y, Z);
|
||||
|
||||
% Customize the plot
|
||||
xlabel('x'); % Label for x-axis
|
||||
ylabel('y'); % Label for y-axis
|
||||
zlabel('f(x, y)'); % Label for z-axis
|
||||
title(plot_title, 'Interpreter', 'latex', 'FontSize', 16); % Title of the plot
|
||||
colorbar;
|
||||
|
||||
% save the figure
|
||||
if strcmp(filename, '') == 0
|
||||
print(gcf, filename, '-dpng', '-r300');
|
||||
end
|
||||
end
|
||||
@@ -0,0 +1,36 @@
|
||||
function plotContour(fun, x_lim, y_lim, size, plot_title, filename)
|
||||
% plot the contour of a function
|
||||
% fun: The function to 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
|
||||
%
|
||||
|
||||
global image_width,
|
||||
global image_height;
|
||||
|
||||
% Generate a grid for x and y
|
||||
x_space = linspace(x_lim(1), x_lim(2), size);
|
||||
y_space = linspace(y_lim(1), y_lim(2), size);
|
||||
[X, Y] = meshgrid(x_space, y_space);
|
||||
|
||||
% Evaluate the function on the grid
|
||||
Z = fun(X, Y);
|
||||
|
||||
% Contour
|
||||
figure('Name', 'Contours of f(x,y)', 'NumberTitle', 'off');
|
||||
set(gcf, 'Position', [100, 100, image_width, image_height]); % Set the figure size
|
||||
contour(X, Y, Z);
|
||||
|
||||
% Customize the plot
|
||||
xlabel('x'); % Label for x-axis
|
||||
ylabel('y'); % Label for y-axis
|
||||
title(plot_title, 'Interpreter', 'latex', 'FontSize', 16); % Title of the plot
|
||||
colorbar;
|
||||
|
||||
% save the figure
|
||||
if strcmp(filename, '') == 0
|
||||
print(gcf, filename, '-dpng', '-r300');
|
||||
end
|
||||
end
|
||||
@@ -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
|
||||
@@ -0,0 +1,28 @@
|
||||
function plotItersOverGamma(gamma, iterations, plot_title, filename)
|
||||
% 3D plots a function
|
||||
% fun: 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;
|
||||
|
||||
figure('Name', 'Iterations_over_gamma', 'NumberTitle', 'off');
|
||||
set(gcf, 'Position', [100, 100, image_width, image_height]); % Set the figure size
|
||||
plot(gamma, iterations, '*r', 'LineWidth', 2);
|
||||
|
||||
% Customize the plot
|
||||
title(plot_title, 'Interpreter', 'latex', 'FontSize', 16); % Title of the plot
|
||||
xlabel('\gamma') ;
|
||||
ylabel('Iterations');
|
||||
|
||||
% save the figure
|
||||
if strcmp(filename, '') == 0
|
||||
print(gcf, filename, '-dpng', '-r300');
|
||||
end
|
||||
end
|
||||
@@ -0,0 +1,43 @@
|
||||
function plotPointsOverContour(points, contour_fun, x_lim, y_lim, size, 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;
|
||||
|
||||
% Generate a grid for x and y
|
||||
x_space = linspace(x_lim(1), x_lim(2), size);
|
||||
y_space = linspace(y_lim(1), y_lim(2), size);
|
||||
[X, Y] = meshgrid(x_space, y_space);
|
||||
|
||||
% Evaluate the function on the grid
|
||||
Z = contour_fun(X, Y);
|
||||
|
||||
% 2D plot
|
||||
figure('Name', '(x,y) convergence', 'NumberTitle', 'off');
|
||||
set(gcf, 'Position', [100, 100, image_width, image_height]); % Set the figure size
|
||||
plot(points(:, 1), points(:, 2), '-or');
|
||||
hold on
|
||||
contour(X, Y, Z);
|
||||
% Customize the plot
|
||||
xlim(x_lim);
|
||||
ylim(y_lim);
|
||||
xlabel('x'); % Label for x-axis
|
||||
ylabel('y'); % Label for y-axis
|
||||
grid on
|
||||
|
||||
title(plot_title, 'Interpreter', 'latex', 'FontSize', 16); % Title of the plot
|
||||
colorbar;
|
||||
|
||||
% save the figure
|
||||
if strcmp(filename, '') == 0
|
||||
print(gcf, filename, '-dpng', '-r300');
|
||||
end
|
||||
end
|
||||
Binary file not shown.
Reference in New Issue
Block a user