\section{Euler-Lagrange Formalism}

The Swiss mathematician Leonhard Euler and the Italo-French mathematician Joseph-Louis Lagrange (see also figure \ref{dithmar_euler-lagrange}) developed the Euler-Lagrange theory in the 1750s. Actually motivated from mechanics, it nowadays pops up in many different fields, e.g.\ field theory or functional minimisation. Also, there exist connections to iterative algorithms like Gradient Descent or Newton's Method.
\begin{figure}[htb]
	\centering
	\includegraphics[width=0.3\textwidth]{dithmar_euler.pdf} \quad
	\includegraphics[width=0.3\textwidth]{dithmar_lagrange.pdf} \\
	\caption{\small \textbf{Left:} Leonhard Euler (1707-1783). \textbf{Right:} Joseph-Louis Lagrange (1736-1813). \textbf{Source:} Wikipedia.}
	\label{dithmar_euler-lagrange}
\end{figure}

\subsection{Basic Observations}
Besides the standard calculus, which we know from school, there also exists the calculus of variations. For both concepts, let's briefly sketch what they have in common and what is different for them.

\textbf{Standard Calculus} \\
In standard calculus, one considers real-valued functions $f: \mathbb{R}^N \to \mathbb{R}; (\xi_1, \dots, \xi_N) \in \mathbb{R}^N$. We know from school that $f$ has a minimum in a point $\xi$ if its first derivative $f'(\xi) = 0$ which is the necessary condition for a minimiser (1-D case). Furthermore, if $f$ is strictly convex and $f'(\xi) = 0$, then $\xi$ is a unique minimiser of $f$.

\textbf{Calculus of Variations} \\
The calculus of variations considers real-valued "functions" that do not take one or several variables as arguments, but that take a function as argument. In this case, one talks about \textit{functionals}. Such a functional $\mathcal{E}: \mathcal{F} \to \mathbb{R}$ is minimised by a function $u$ if $u$ satisfies the so-called \textit{Euler-Lagrange equations}. This is -- like the first derivative in standard calculus -- a necessary condition for a minimiser. Similarly, if $\mathcal{E}$ is strictly convex and satisfies the Euler-Lagrange equations, then $u$ is a unique minimiser of $\mathcal{E}$.

\textbf{Functionals} \\
For the rest of this paragraph, we want to consider functionals that have the following structure
\begin{equation}
	\label{dithmar_energy-functional}
	\mathcal{E}(u) = \int_a^b \underbrace{(u-f)^2}_{\mbox{data term}} + \alpha \underbrace{\Psi_S(|\nabla u|^2)}_{\mbox{smoothness term}} \; \mathrm d\textbf{x}
\end{equation}
where $u$ is the (unknown) function for which we minimise, $\alpha$ is a smoothness weight which controls the influence of the smoothness term and $\Psi_S$ is a regulariser that penalises deviations from smoothness. Please note that one can also use a penaliser function $\Psi_D$ to robustify the data term against outliers. The concept of penalising deviations from (measured) data was also part of "Approximation and Fitting I" that was presented by Laurent Hoeltgen.

\subsection{Euler-Lagrange equations}
Let's consider the 1-D and the 2-D case of the Euler-Lagrange theorem.

\subsubsection{Theorem (1-D case)}
A \textit{smooth function} $u(x)$ with $x \in [a,b]$ that minimises the 1-D energy functional
\begin{equation}
	\mathcal{E}(u) = \int_a^b \mathcal{L}(x,u,u') \; \mathrm dx
\end{equation}
satisfies necessarily the Euler-Lagrange equation
\begin{equation}
	\label{dithmar_EL1d}
	\mathcal{L}_u - \frac{\mathrm d}{\mathrm dx} \mathcal{L}_{u'} = 0
\end{equation}
and the so-called natural boundary conditions
\begin{equation}
	\mathcal{L}_{u'} = 0
\end{equation}
for the boundary ($x = a$ and $x = b$) where $\mathcal{L}$ is sometimes called the Lagrangian. \\
Informally speaking, \textit{smooth function} means that one "doesn't get into trouble" when making computations, or more formally, the function $u$ should be as often differentiable as needed.

\subsubsection{Theorem (2-D case)}
The 2-D case is a simple extension of the 1-D scenario. In principle, $u$ takes two variables as arguments, $x_1$ and $x_2$. This requires to use partial derivatives instead of "normal" ones. The theorem reads as follows: \\
A minimiser of the 2-D energy functional
\begin{equation}
	\mathcal{E}(u) = \int_\Omega \mathcal{L}(x_1,x_2,u,u_{x_1},u_{x_2}) \; \mathrm d\textbf{x}
\end{equation}
satisfies necessarily the Euler-Lagrange equation
\begin{equation}
	\mathcal{L}_u - \frac{\partial}{\partial x_1} \mathcal{L}_{u_{x_1}} - \frac{\partial}{\partial x_2} \mathcal{L}_{u_{x_2}} = 0
\end{equation}
with the natural boundary conditions
\begin{equation}
	\textbf{n}^\top \left(\begin{array}{c}\mathcal{L}_{u_{x_1}} \\ \mathcal{L}_{u_{x_2}}\end{array}\right) = 0
\end{equation}
at the image boundaries $\partial\Omega$ with the normal vector $\textbf{n}$.

\subsubsection{Example for a Variational Problem}
Let's give a short example that illustrates how the Euler-Lagrange equations can be used. The task is to find the shortest plane curve joining two points $A$ and $B$. One can reformulate this as a variational problem as follows: \\
Find the curve $u(x)$ for which the functional
\begin{equation}
	\mathcal{E}(u(x)) = \int_a^b \sqrt{1 + \left(\frac{\mathrm d}{\mathrm dx} u(x)\right)^2} \; \mathrm dx
\end{equation}
achieves its minimum where $a \leq x \leq b$. \\
We know that the solution to this problem is a straight line segment joining $A$ and $B$. The explanation can be given using the Euler-Lagrange formalism. We have the Lagrangian
\begin{equation}
	\mathcal{L} = \sqrt{1 + \left(\frac{\mathrm d}{\mathrm dx} u(x)\right)^2},
\end{equation}
and we need the following ingredients where the more convenient notation $u'$ is used to indicate the derivative of $u$ with respect to $x$:
\begin{eqnarray}
	\mathcal{L}_u & = & 0 \\
	\nonumber \mathcal{L}_{u'} & = & \frac{1}{2 \cdot \sqrt{1 + (u')^2}} \cdot 2 \cdot u' \\
	& = & \frac{u'}{\sqrt{1 + (u')^2}}
\end{eqnarray}
We can now easily see that $u'$ has to be a constant to satisfy the Euler-Lagrange equation
\begin{equation}
	\frac{\mathrm d}{\mathrm dx}\frac{u'}{\sqrt{1 + (u')^2}} = 0.
\end{equation}
$u'$ being a constant means in other words that $u(x)$ has to be a straight line.

\subsection{First and Second Variation of an Energy Functional}
So far, we assumed that we are able to compute the derivatives $\mathcal{L}_u$, $\mathcal{L}_{u_{x_1}}$ and $\mathcal{L}_{u_{x_2}}$. The problem in this case is that $u$ and its derivatives are functions, not values. But how do we compute a derivative with respect to a function? To answer this question, let's first introduce some notation where we distinguish the continuous and the discrete setting. In the continuous setting, we refer to a functional by using $\mathcal{E}$. The functional $\mathcal{E}$ takes functions $u$ as arguments. In the discrete setting, instead, we have a function $E$ that takes several variables -- represented by a vector $\textbf{u}$ -- as arguments. \\
Now we have the basic ingredients to answer the question how to compute a derivative with respect to a function. The goal is to construct "something" of which we can compute the derivative. This is achieved by actually two steps: a function that brings us to the function space $\mathcal{F}$ and a function that brings us back to $\mathbb{R}$:
\begin{itemize}
	\item $\mathcal{S}: \mathbb{R} \to \mathcal{F}; \; \varepsilon \mapsto u + \varepsilon v$ brings us to the function space $\mathcal{F}$.
	\item $\mathcal{E}: \mathcal{F} \to \mathbb{R}$ brings us back to $\mathbb{R}$.
	\item[$\Rightarrow$] $\mathcal{P}: \mathbb{R} \to \mathbb{R}; \; \varepsilon \mapsto \mathcal{E} (u + \varepsilon v)$.
\end{itemize}
Let's think of $u$ and $v$ being vectors $\textbf{u}$ and $\textbf{v}$, as we do for the discrete setting. In this case we can consider the set of points
\begin{equation}
	\{\textbf{u} + \varepsilon \textbf{v} | \varepsilon \in \mathbb{R}\}
\end{equation}
which describes a straight line along the vector $\textbf{v}$ through $\textbf{u}$. The derivative of
\begin{equation}
	\varepsilon \mapsto E(\textbf{u} + \varepsilon \textbf{v})
\end{equation}
at $\varepsilon = 0$ is described by
\begin{equation}
	\label{dithmar_discrete-1st-directional-derivative}
	\frac{\partial}{\partial \varepsilon} E(\textbf{u} + \varepsilon \textbf{v})\Big|_{\varepsilon = 0}.
\end{equation}
The latter is called directional derivative. Equation \eqref{dithmar_discrete-1st-directional-derivative} and
\begin{equation}
	\frac{\partial^2}{\partial \varepsilon^2} E(\textbf{u} + \varepsilon \textbf{v})\Big|_{\varepsilon = 0}
\end{equation}
can also be formulated in the continuous setting, resulting in the first and second variation, respectively:
\begin{eqnarray}
	\frac{\partial}{\partial \varepsilon} \mathcal{E}(u + \varepsilon v)\Big|_{\varepsilon = 0}, \\
	\frac{\partial^2}{\partial \varepsilon^2} \mathcal{E}(u + \varepsilon v)\Big|_{\varepsilon = 0}.
\end{eqnarray}
The intuition should tell us that the first variation is something like a first derivative and the second variation something like a second derivative. Since we are heading for iterative algorithms, we consider from now on only the discrete case. \\
This discrete version of the first variation is given by
\begin{eqnarray}
	\nonumber \textbf{D}E(\textbf{u})\textbf{v} & = & \frac{\partial}{\partial \varepsilon} E(\textbf{u} + \varepsilon \textbf{v})\Big|_{\varepsilon = 0} \\
	& = & \nabla E(\textbf{u}) \cdot \textbf{v},
\end{eqnarray}
and of the second variation by
\begin{eqnarray}
	\nonumber \textbf{D}^2E(\textbf{u})\textbf{v} & = & \left. \frac{\partial^2}{\partial \varepsilon^2} E(\textbf{u} + \varepsilon \textbf{v}) \right|_{\varepsilon = 0} \\
	\nonumber & = & \left. \frac{\partial}{\partial \varepsilon} (\nabla E(\textbf{u} + \varepsilon \textbf{v}) \cdot \textbf{v}) \right|_{\varepsilon = 0} \\
	\nonumber & = & \left. \frac{\partial}{\partial \varepsilon} \left(\left(\frac{\partial}{\partial u_1} E(\textbf{u} + \varepsilon \textbf{v}), ..., \frac{\partial}{\partial u_N} E(\textbf{u} + \varepsilon \textbf{v})\right) \cdot \textbf{v}\right) \right|_{\varepsilon = 0} \\
	\nonumber & = & \left(\left(\left(\frac{\partial^2}{\partial u_1 \partial u_i} E(\textbf{u})\right), ..., \left(\frac{\partial^2}{\partial u_N \partial u_i} E(\textbf{u})\right)\right)_{i=1,...,N} \cdot \textbf{v}\right) \cdot \textbf{v} \\
	\nonumber & = & \left(\left(\underbrace{\left(\frac{\partial^2}{\partial u_j \partial u_i} E(\textbf{u})\right)}_{\mathcal{H} E(\textbf{u})}\right)_{i,j=1,..,N} \cdot \textbf{v}\right) \cdot \textbf{v} \\
	& = & \textbf{v}^\top \mathcal{H} E(\textbf{u}) \textbf{v}.
\end{eqnarray}
This can now be used in iterative algorithms.

\subsection{Appearance of the Euler-Lagrange Formalism}
In this part we are interested in iterative algorithms where the Euler-Lagrange formalism comes into play. For this purpose, we recall
\begin{eqnarray}
	\nabla E(\textbf{u}) & = & \left(\frac{\partial}{\partial u_j}(E(\textbf{u}))\right)_{j=1,\dots,N}, \\
	\mathcal{H} E(\textbf{u}) & = & \left(\frac{\partial^2}{\partial u_j \partial u_i}(E(\textbf{u}))\right)_{i,j=1,\dots,N}.
\end{eqnarray}

\subsubsection{Gradient Descent Algorithm}
The intuition for the Gradient Descent Algorithm is simple: $\nabla E(\textbf{u})$ itself points always in the direction of the steepest ascent where $-\nabla E(\textbf{u})$ obviously points in the opposite direction, i.e.\ the steepest descent. Since we want to do a minimisation, $-\nabla E(\textbf{u})$ gives us the search direction for each iteration step of the algorithm:
\\ \ \\
\verb!give a starting point! $x$ \\
\verb!repeat! \\
\verb!  ! $\textbf{d} = -\nabla E(\textbf{u})$ \\
\verb!  Choose step size! $\gamma$ \verb!(e.g. by backtracking line search)! \\
\verb!  ! $\textbf{x} := \textbf{x} + \gamma \cdot \textbf{d}$ \\
\verb!until stopping criterion is satisfied!
\\ \ \\
The algorithm can also be illustrated as shown in figure \ref{dithmar_gradient-descent}.
\begin{figure}[htb]
	\centering
	\includegraphics[width=0.35\textwidth]{dithmar_gradient-descent.png} \\
	\caption{\small Illustration of gradient descent. \textbf{Source:} Wikipedia.}
	\label{dithmar_gradient-descent}
\end{figure} 

\subsubsection{Newton's Method}
We also want to briefly have a look at Newton's Method:
\\ \ \\
\verb!give a starting point! $x$ \\
\verb!repeat! \\
\verb!  ! $\textbf{d} = -[\mathcal{H}E(\textbf{u})]^{-1} \nabla E(\textbf{u})$ \\
\verb!  ! $\textbf{x} := \textbf{x} + \textbf{d}$ \\
\verb!until stopping criterion is satisfied!
\\ \ \\
It is important to mention that in a optimisation problem a strictly convex function is desirable. This notion is equivalent to a positive definite Hessian matrix. But what can we do if the Hessian is not positive definite? A solution to this problem could be to perturb the Hessian until it is positive definite which is done in the Trust-Region Method.

\subsection{Connection to the Euler-Lagrange Formalism}
We will now have a look at the connection between the iterative algorithms briefly mentioned above and the Euler-Lagrange formalism. For this purpose we assume the 1-D version of \eqref{dithmar_energy-functional}, i.e.
\begin{equation}
	\label{dithmar_energy-functional-connection}
	\mathcal{E}(u) = \int_a^b (u-f)^2 + \alpha \Psi_S((u')^2) \; \mathrm dx
\end{equation}
This functional is obviously continuous. The question is now how to get a discrete version for the usage in iterative algorithms out of it. Two approaches are possible:
\begin{itemize}
	\item First employ the necessary optimality condition, then discretise. (OD)
	\item First discretise, then employ the necessary optimality condition. (DO)
\end{itemize}

\subsubsection{Necessary Optimality Condition $\to$ Discretise: OD}
First employing the necessary optimality conditions, i.e.\ employing the Euler-Lagrange equation, gives
\begin{eqnarray}
	\nonumber 0 & = & 2 (u-f) - \frac{\mathrm d}{\mathrm dx} \Big[2 \alpha \Psi_S'\big((u')^2\big) \cdot u'\Big] \\
	& = & u-f - \alpha \frac{\mathrm d}{\mathrm dx} \Big[\Psi_S'\big((u')^2\big) \cdot u'\Big]
\end{eqnarray}
After discretisation, we end up in
\begin{eqnarray}
	\nonumber 0 & = & u_j - f_j - \alpha \frac{\mathrm d}{\mathrm dx} \Big[\Psi_S'\big((u')^2\big) \cdot u'\Big]\Big|_{x=x_j} \\
	\nonumber & = & u_j - f_j - \alpha \Psi_S'\left(\frac{(u_{j+1} - u_j)^2}{\Delta x^2}\right) \cdot \frac{u_{j+1} - u_j}{\Delta x^2} \\
	\label{dithmar_od} & & \quad \quad \quad \; \: + \alpha \Psi_S'\left(\frac{(u_j - u_{j-1})^2}{\Delta x^2}\right) \cdot \frac{u_j - u_{j-1}}{\Delta x^2} \quad , j = 1,...,N.
\end{eqnarray}

\subsubsection{Discretise $\to$ Necessary Optimality Condition: DO}
First discretising the energy functional gives
\begin{equation}
	E(\textbf{u}) = \sum_{j=1}^N (u_j - f_j)^2 + \alpha \Psi_S\Big((\textbf{u}')^2\Big)\Big|_{x=x_j}
\end{equation}
After that, wen can employ the necessary conditions for a minimiser, i.e.\ all partial derivatives must be 0:
\begin{eqnarray}
	\nonumber 0 & = & \frac{\partial}{\partial u_j} \left(E(\textbf{u})\right) \\
	\nonumber & = & u_j - f_j - \alpha \Psi_S'\left(\frac{(u_{j+1} - u_j)^2}{\Delta x^2}\right) \cdot \frac{u_{j+1} - u_j}{\Delta x^2} \\
	\label{dithmar_do} & & \quad \quad \quad \; \: + \alpha \Psi_S'\left(\frac{(u_j - u_{j-1})^2}{\Delta x^2}\right) \cdot \frac{u_j - u_{j-1}}{\Delta x^2} \quad , j = 1,...,N.
\end{eqnarray}
As one can easily verify, \eqref{dithmar_od} and \eqref{dithmar_do} are exactly the same.

\subsection{Examples in Image Processing}
In the following, we want to have a look at two different examples from image processing: removing Gaussian noise and removing salt-and-pepper noise.

\subsubsection{Removing Gaussian Noise}
For this example, let's consider again the energy functional
\begin{equation}
	\mathcal{E}(u) = \int_a^b (u-f)^2 + \alpha \Psi_S(|\nabla u|^2) \; \mathrm d\textbf{x}
	\label{dithmar_example-gaussian-noise-energy-functional}
\end{equation}
with the Charbonnier regulariser
\begin{equation}
	\Psi_S(|\nabla u|^2) = 2(\lambda^2 \sqrt{1 + |\nabla u|^2 / \lambda^2} - \lambda^2)
\end{equation}
For the further computation, OD was applied to discretise \eqref{dithmar_example-gaussian-noise-energy-functional} for the usage in an iterative algorithm with explicit time stepping where further details are omitted. Let's compare the results after 20 iterations which are given in figure \ref{dithmar_example1}.
\begin{figure}
	\begin{center}
		\includegraphics[width=0.3\textwidth]{dithmar_couple.pdf} \quad
		\includegraphics[width=0.3\textwidth]{dithmar_couple40.pdf} \\[0.2cm]
		\includegraphics[width=0.3\textwidth]{dithmar_couple40_1_20.pdf} \quad
		\includegraphics[width=0.3\textwidth]{dithmar_couple40_2_20.pdf} \quad
		\includegraphics[width=0.3\textwidth]{dithmar_couple40_5_20.pdf}
	\end{center}
	\caption{\small \textbf{Top Left:} Original image. \textbf{Top Right:} Noisy version with Gaussian noise ($\sigma = 40$). \textbf{Bottom row:} Denoised version with $\lambda = 1,2,5$. The higher the contrast parameter $\lambda$, the more noise is removed, but the image is also more blurred.}
	\label{dithmar_example1}
\end{figure}

\subsubsection{Removing Salt-and-Pepper Noise}
A second example is the removal of salt-and-pepper noise given in figure \ref{dithmar_example2}. This example is taken from the DIC lecture.
\begin{figure}
	\begin{center}
		\includegraphics[width=0.35\textwidth]{dithmar_sp_noisy.pdf} \quad \quad
		\includegraphics[width=0.35\textwidth]{dithmar_sp_noisefree.pdf}
	\end{center}
	\caption{\small \textbf{Left:} Noisy version with salt-and-pepper noise. \textbf{Right:} Denoised version using the regularised data term $\Psi_D(s^2) = 2 (\sqrt{s^2 + \varepsilon^2} - \varepsilon)$ and a Charbonnier regulariser as smoothness term $\Psi_S(s^2) = 2 (\lambda^2 \sqrt{1 + s^2 / \lambda^2} - \lambda^2)$. \textbf{Author:} J. Weickert.}
	\label{dithmar_example2}
\end{figure}

\subsection{Summary}
We have seen that the Euler-Lagrange equations are a necessary condition for a minimiser of an energy functional. They are used in several approaches like iterative algorithms. There we have seen the Gradient Descent Algorithm and Newton's Method which are suitable to find a local extremum. For those algorithms, the Euler-Lagrange equations provide the search direction $\textbf{d}$ (cf. figure \ref{dithmar_E-L}).
\begin{figure}[htb]
	\centering
	\includegraphics[width=0.45\textwidth]{dithmar_E-L.pdf}
	\caption{\small The Euler-Lagrange equations provide the search direction for either the Gradient Descent Algorithm or Newton's Method.}
	\label{dithmar_E-L}
\end{figure}
Furthermore, we have seen that it doesn't matter whether we employ first the necessary optimality condition -- i.e.\ the Euler-Lagrange equations -- and then discretise or vice versa. Examples showed us that energy functional minimisation is a possibility to perform image denoising.

\subsection{References}
\begin{itemize}
	\item Pablo Pedregal. Introduction to Optimization. Springer New York. 2004.
	\item I.M. Gelfand, S.V. Fomin. Calculus of Variations. Dover New York. 2000.
	\item Euler-Lagrange equations (Wikipedia). \\
		\texttt{http://en.wikipedia.org/wiki/Euler-Lagrange\_equation}
	\item Lecture "Differential Equations in Image Processing and Computer Vision". SS 2008.
	\item Lecture "Numerical Algorithms for Visual Computing III". SS 2009.
	\item L. Bar, G. Sapiro. Generalized Newton-Type Methods for Energy Formulations in Image Processing. SIAM Journal on Imaging Sciences, 2009. \\
		\texttt{http://www.ece.umn.edu/users/barxx002/papers/Siam09.pdf}
\end{itemize}