\chapter{Background}

This section provides the necessary background for understanding the theoretical aspects of this project. We begin by introducing the notion of exact computation \cite{ref3}, a new emerging paradigm in the field of modern computation. We also discuss the significance of this work in the context of exact geometric computation \cite{ref4}. Finally, we introduce the main concepts of optimisation and linear programming that the reader should be familiar with.

\section{Exact Computation}

The underlying nature of computation is par excellence numerical: numbers have been at the heart of all calculations since very ancient times, while modern mathematics developed a whole range of theories to formalise many aspects of computable numbers. Early computers had the sole purpose of performing large complex calculations (the so-called number crunchers), which then turned into the original mass-produced computers, in the form of friendly pocket calculators. Although computers have evolved significantly and moved towards more abstracted and higher-level programs, numerical computation remains a major pillar of modern computer technology.
\\

In particular, scientific computation is a rapidly growing inter-disciplinary field that makes use of advanced numerical capabilities. It is considered as adding a new dimension to the classical methods of theory and experimentation, sometimes referred to as the ``third scientific method'' of computation \cite{ref3}.

\subsection{Fixed-Point Paradigm}

At its core, scientific computation is subject to the \textit{fixed-precision paradigm} of computation. Under this representation, numbers are expressed using a fixed number of digits after the decimal/radix point, thus providing fixed computational precision (usually machine-dependant). 
\\

Based on this specification, one can use several approaches to limit the unavoidable rounding errors. For instance, a \textit{mild form} of fixed-precision can be applied such that computations are performed up to a user-specified precision level. Although we can set a very high level of precision that we may think is satisfactory, the build-up of \textit{round-off errors} that accumulate may produce totally unexpected results. We can illustrate this with an example from \cite{ref5}. Consider the sequence $a_n$ defined recursively as:

\[
a_n = \left\{\begin{array}{cll}
\displaystyle\frac {11}2    & , & n = 0 \\\\
\displaystyle\frac {61}{11} & , & n = 1 \\\\
111 - \displaystyle\frac{1130 - \frac {3000}{a_{n-1}}}{a_n} & , & n \ge 2
\end{array}\right. .
\]
\\

Using the Unix utility \lit*{bc}, we can compute the terms of the sequence $a_n$ up to some fixed precision of $k$ decimal places, which we shall denote by $a_n^{(k)}$. Performing calculations with 5 decimal places, gives (note that the results have been rounded for presentation purposes, as in \cite{ref5}):
\\

\[\begin{array}{c}
\begin{tabular}{|c|r|}
\hline
$a_0^{(5)}$ & $5.500$  \\ \hline
$a_1^{(5)}$ & $5.500$  \\ \hline
$a_2^{(5)}$ & $5.500$  \\ \hline
$a_3^{(5)}$ & $5.500$  \\ \hline
$a_4^{(5)}$ & $ 5.648$ \\ \hline
$a_5^{(5)}$ & $ 5.242$ \\ \hline
\end{tabular}

\hspace{2cm}

\begin{tabular}{|c|r|}
\hline
$a_6^{(5)}$    & $ -3.241$ \\ \hline
$a_7^{(5)}$    & $  283.1$ \\ \hline
$a_8^{(5)}$    & $103.738$ \\ \hline
$a_9^{(5)}$    & $100.209$ \\ \hline
$a_{10}^{(5)}$ & $100.012$ \\ \hline
$a_{11}^{(5)}$ & $100.001$ \\ \hline
\end{tabular}
\end{array}\]
\\

One would therefore believe that the sequence converges to $100$. However, computing the number $a_{100}$ with higher and higher precisions will contradict our expectations (the ``exponents'' below indicate repeating digits, e.g. $1.2^34=1.2224$):

\[\begin{array}{c}
\begin{tabular}{|c|l|}
\hline
$a_{100}^{(5)}$   & $100.0^41$           \\ \hline
$a_{100}^{(30)}$  & $100.0^{29}1$        \\ \hline
$a_{100}^{(60)}$  & $100.0^{57}997$      \\ \hline
$a_{100}^{(100)}$ & $100.0^{17}98\ldots$ \\ \hline
\end{tabular}

\hspace{2cm}

\begin{tabular}{|c|l|}
\hline
$a_{100}^{(110)}$ & $100.0^{7}92$        \\ \hline
$a_{100}^{(120)}$ & $-3.790\ldots$       \\ \hline
$a_{100}^{(130)}$ & $5.9^{7}8697\ldots$  \\ \hline
$a_{100}^{(140)}$ & $5.9^{7}87925\ldots$ \\ \hline
\end{tabular}
\end{array}\]
\\

Here we can notice that if the precision is either $5$, $30$, $60$, $100$ or even $110$ decimal places, then our expectation from above holds ($a_n\to 100$ as $n$ grows larger). However, when increasing the precision even further, we obtain rather spurious results. 

The actual limit of the sequence is equal to $6$, since we can evaluate the general term $a_n$ as:

\[a_n = \frac {6^{n+1} + 5^{n+1}}{6^n + 5^n}.\]
\\

This example therefore shows how even the \textit{mild form} of fixed-precision can lead to flawed results. It is therefore unclear what level of precision needs to be set in advance to a program such as \lit*{bc} in order to obtain the correct answer. As we could see, up until 110 decimal places we got roughly the same (wrong) result and we actually had to consider precision above 130 decimal places to arrive at the right answer.

\subsection{Floating-Point Representation}

Similar to the fixed precision paradigm, modern computers use floating-point arithmetic to perform real number calculations. This representation can be seen as a trade-off between range and precision when approximating real number: with a given precision, the floating-point model is able to represent both numbers of small magnitude with many bits of significance or conversely, large magnitude and few bits of significance. While many different representations have been developed in the past, the one defined by IEEE 754 has emerged as the industry standard. This standard is a step forward towards addressing the issue of portability and therefore makes errors in floating point computation \textit{machine-independent}.
\\

However, this representation is essentially just as inaccurate as the fixed-point model, since we must approximate real numbers by their nearest representable one. Thus, the same rounding errors will occur in practice, e.g. when small inaccuracies propagate in successive iterations like the one that we have just seen when computing a simple mathematical limit.

\subsection{Towards an Alternative to the f.p. Paradigm} 

Even with industry-leading standards, rather intractable problems arise from the presence of round-off errors and compromise the \textit{robustness} of the f.p. paradigm. We can, at the arithmetic level, increase precision via techniques such as double extended precision, guard-bits, gradual underflow etc, which can usually be implemented in hardware. Interval arithmetic is another well-known method. Looking from a geometric perspective, an idea would be to divide the input and computed data into combinatorial and numerical, and to give precedence to the former when making decisions. An argument in favour of this approach is that we can allow the numerical data to be perturbed in order to maintain the combinatorial data. This avoids ``topological inconsistencies'' and can be implemented for simple cases, but the intractable nature of combinatorial problems makes it rather difficult to deal with more general cases \cite{ref3}.
\\

Therefore, in the light of these non-robustness issues, one needs a completely different approach to handle cases in which the correct and exact answer is required. A new direction in the literature of computation is called, unsurprisingly, the \textit{exact computation paradigm}. According, to \cite{ref3} or \cite{ref4}, this paradigm assumes a computational process that:

\begin{enumerate}
	\item represents all the underlying mathematical objects \textit{exactly};
	\item all branching decisions are error-free.
\end{enumerate}

As a result, multi-precision arithmetic is a necessary condition (but not sufficient) for exact computation. A different issue is that exact computation will naturally come at the expense of performance. It therefore makes sense to target only those applications which are not \textit{cycle-critical} (i.e. we afford to incur some sort of computational slow-down). For example, exact computation cannot be avoided in computational number theory and in many aspects of algebra (e.g. testing the irreducibility of a polynomial).
\\

Finally, we conclude this section with the idea of \textit{weak exact computation}, similar in intent with the mild form of the f.p. paradigm. As explained in the numerical example above, one starts the computations using some fixed bound $k$ on the precision. However, one would need to perform a thorough analysis as to what minimal values of $k$ will indeed give the exact results. This is clearly a non-trivial problem and several suggestions have been proposed: in the same \cite{ref3}, the theory of root bounds has been developed; in \cite{ref5}, the exact real arithmetic offers lower and upper bounds guarantees that are trustworthy.

\section{Linear Programming}

A great deal of problems encountered in the real world involve maximisation or minimisation of certain quantities. Most often, these quantities we seek to optimise are profits (in the case of maximisation) as well as costs (in the case of minimisation). In linear programming we therefore aim to minimise/maximise a certain linear function (which we usually call objective function) subject to some linear constraints that describe the restrictions of our problem. 
\\

One would immediately think that calculus, developed by Leibniz and Newton in the 17th century, can deal with this types of problems very elegantly given the arsenal of techniques readily available. However, calculus is inadequate since it can only imply that the maxima and/or minima of some objective function lie on the boundaries of the sets determined by the constraints. Thus, a special set of techniques and algorithms is needed to deal with such linear programming problems (this is what the term ``programming'' really refers to in this case).
\\

As such, linear programming plays a very important role in the fields of mathematics and business, finding many applications in management science and operations research.

\subsection{Linear Programming in Standard Form}

A linear program (LP) is an optimisation problem in which the \textit{objective function} (i.e. the quantity we want to minimise/maximise) is linear in the unknowns and the constraints consist of linear equalities and/or linear inequalities. The exact form of these constraints may be problem-dependant, but one can always rewrite a linear program in the following \textit{standard form}:

\begin{equation} \label{standard-form}
\begin{array}{llllllllll}
\text{Minimise}   & c_1x_1    & + & c_2x_2    & + & \ldots & + & c_nx_n    &     &       \\
\text{subject to} & a_{11}x_1 & + & a_{12}x_2 & + & \ldots & + & a_{1n}x_n & = & b_1   \\
	              & a_{21}x_1 & + & a_{22}x_2 & + & \ldots & + & a_{2n}x_n & = & b_2   \\
	              & & & & \ \vdots & & & & & \\
	              & a_{m1}x_1 & + & a_{m2}x_2 & + & \ldots & + & a_{mn}x_n & = & b_m \\
	              & \mathrlap{x_1, x_2, \ldots x_n \ge 0}. & & & & & & & &   
\end{array}
\end{equation}

\noindent where the coefficients $b_i, c_j, a_{ij}$, $i=\overline{1, m}$, $j = \overline{1, n}$ are fixed real constants and the \textit{decision variables} $x_i$, $i = \overline{1, n}$, are yet to be found. For simplicity, we will assume that all $b_i\ge 0$ (if necessary, one can multiply by $-1$ the equations for which $b_i \le 0$). Lastly, the first $m$ constraints are said to be \textit{main constraints}, while the second $n$ constraints are called \textit{non-negativity constraints}.
\\

We can rewrite the standard form given by \ref{standard-form} using matrix notation to derive the following compact statement:

\begin{equation} \label{standard-form-compact}
\begin{array}{ll}
\text{Minimise}   & \bm{c}^\top \bm{x} \\
\text{subject to} & \bm{Ax} = \bm{b} \\
	              & \bm{x} \ge \bm{0},
\end{array}\\[0.2cm]
\end{equation}

\noindent where $\bm{x}$ is an $n$-dimensional column vector, $\bm{c}^\top$ is an $n$-dimensional row vector, $\bm{A}$ is an $m\times n$ matrix and $\bm{b}$ is an $m$-dimensional column vector with $\bm{b} \ge \bm{0}$. Inequalities of the type $\bm{x}\ge \bm{0}$ are understood to hold component-wise.
\\

To illustrate the point that we can convert any LP problem into standard form, we need to consider various scenarios that may arise in practice. Apart from the fact that we may have equality constraints with negative right hand sides (which we can mitigate upon multiplication with $-1$, as previously mentioned), general LP problems can: 

\begin{enumerate}
	\item be maximisation (instead of minimisation) problems; 
	\item have inequality (instead of equality) constraints;
	\item have free (instead of non-negative) decision variables.
\end{enumerate}

In order to deal with maximisation problems, we need to observe that we can simply invert the objective function and then formulate the corresponding LP problem as a minimisation problem. We thus state the following: 

\begin{theorem} \label{min-max}
Let $f:\Omega \to\mathbb{R}$ be a real valued function and assume that both the minimum and maximum of $f$ are attained within the set $\Omega$. Then:
\[\min_{\bm{x}\in\Omega} f(\bm{x}) = - \max_{\bm{x}\in\Omega} -f(\bm{x}).\\[0.2cm]\]	
\end{theorem}

\begin{proof}
The minimum of $f$ satisfies:

\[\min f(\bm{x}) = f^\star \le f(\bm{x}), \forall\ x\in\Omega.\\[0.2cm]\]

\noindent Similarly, the maximum of $-f$ satisfies:

\[\max -f(\bm{x}) = F^\star \ge -f(\bm{x}), \forall\ x\in\Omega.\\[0.2cm]\]

\noindent But if $-F^\star < f^\star$ then there exists an $\bm{\overline{x}}\in\Omega$ such that $f(\bm{\overline{x}}) < f^\star$ contradicting the optimality of $f^\star$. Conversely, if $-f^\star > F^\star$ then there exists an $\bm{\overline{x}}\in\Omega$ such that $-f(\bm{\overline{x}}) > f^\star$ contradicting the optimality of $F^\star$. In conclusion $f^\star = -F^\star$.
\end{proof}

Now, to address the issue of LP problems with inequality constraints, we need to consider the following problem in which the constraints are given by linear inequalities:

\[
\begin{array}{llllllllll}
\text{Minimise}   & c_1x_1    & + & c_2x_2    & + & \ldots & + & c_nx_n    &     &       \\
\text{subject to} & a_{11}x_1 & + & a_{12}x_2 & + & \ldots & + & a_{1n}x_n & \le & b_1   \\
	              & a_{21}x_1 & + & a_{22}x_2 & + & \ldots & + & a_{2n}x_n & \le & b_2   \\
	              & & & & \ \vdots & & & & & \\
	              & a_{m1}x_1 & + & a_{m2}x_2 & + & \ldots & + & a_{mn}x_n & \le & b_m \\
	              & \mathrlap{x_1, x_2, \ldots x_n \ge 0.} & & & & & & & &   
\end{array}\\[0.2cm]
\]

In this case, one can introduce new positive variables $s_i\ge 0$ (called \textit{slack variables}), where $i = \overline{1, m}$, to convert each inequality constraint into an equality in the following manner:

\[
\begin{array}{llllllllllll}
\text{Minimise}   & c_1x_1    & + & c_2x_2    & + & \ldots & + & c_nx_n    &   &     &     &     \\
\text{subject to} & a_{11}x_1 & + & a_{12}x_2 & + & \ldots & + & a_{1n}x_n & + & s_1 & \le & b_1 \\
	              & a_{21}x_1 & + & a_{22}x_2 & + & \ldots & + & a_{2n}x_n & + & s_2 & \le & b_2 \\
	              & & & & \ \vdots & & & & & \\
	              & a_{m1}x_1 & + & a_{m2}x_2 & + & \ldots & + & a_{mn}x_n & + & s_m & \le & b_m \\
	              & \mathrlap{x_1, x_2, \ldots x_n \ge 0,} & & & & & & & & & &                   \\
	              & \mathrlap{s_1, s_2, \ldots s_m \ge 0.} & & & & & & & & & &
\end{array} \\[0.2cm]
\]

It is now easy to see that this is an LP problem in standard form with $m + n$ decision variables and coefficient matrix having the form $[\bm{A}, \bm{I}]$ (where $\bm{A}$ is the original matrix as defined in the standard form \ref{standard-form}, and $\bm{I}$ is the $m\times m$ identity matrix -- corresponding to the slack variables $s_i$ that have been added). In matrix notation, we can write that:

\[
\begin{array}{ll}
\text{Minimise}   & \bm{c}^\top \bm{x} \\
\text{subject to} & \bm{Ax} \le \bm{b} \\
	              & \bm{x} \ge \bm{0},
\end{array} \iff 
\begin{array}{ll}
\text{Minimise}   & \bm{c}^\top \bm{x} \\
\text{subject to} & \bm{Ax} + \bm{s} = \bm{b} \\
	              & \bm{x} \ge \bm{0}, \bm{s} \ge \bm{0},
\end{array}\\[0.2cm]
\]

\noindent where $\bm{s} = [s_1, s_2, \ldots, s_m]^\top$. Hence, the slack variables characterise the difference $\bm{b} - \bm{Ax}$. Similarly, one can deal with reversed inequality constraints such as:

\[a_{i1}x_1 + a_{i2}x_2 + \ldots + a_{in}x_n \ge b_i\\[0.2cm]\]

\noindent by introducing \textit{surplus variables} $s_i\ge 0$ (also known as \textit{excess variables}):

\[a_{i1}x_1 + a_{i2}x_2 + \ldots + a_{in}x_n - s_i = b_i.\\[0.2cm]\]

\noindent Again, the resulting LP problem will have $n + m$ decision variables and we can write it in compact matrix form as follows:

\[
\begin{array}{ll}
\text{Minimise}   & \bm{c}^\top \bm{x} \\
\text{subject to} & \bm{Ax} \ge \bm{b} \\
	              & \bm{x} \ge \bm{0}
\end{array} \iff 
\begin{array}{ll}
\text{Minimise}   & \bm{c}^\top \bm{x} \\
\text{subject to} & \bm{Ax} - \bm{s} = \bm{b} \\
	              & \bm{x} \ge \bm{0}, \bm{s} \ge \bm{0},
\end{array}\\[0.2cm]
\]

\noindent where surplus variables characterise the difference $\bm{Ax} - \bm{b}$.
\\

Finally, let us deal with LP problems in which decision variables $x_i$ are free, i.e. the unknowns can either be positive or negative. There are two different techniques that allow one to convert such a problem into standard form:

\begin{enumerate}[(a)]
	\item If $x_i$ is unconstrained, then we can use the substitution: $x_i = x_i^+ - x_i^-$, where $x_i^+, x_i^-\ge 0$. This means that the LP problem is now in standard form and has $n + 1$ decision variables, namely: $x_1, x_2, \ldots, x_{i-1}, x_i^+$, $x_i^{-}, x_{i+1}, \ldots, x_n$.
	\item If $x_i$ is unconstrained, then we can use any equality constraint to eliminate $x_i$ so that the standardised LP problem will now have $n - 1$ decision variables and $m - 1$ constraints. The value of the free variable $x_i$ can then be determined from the equation used to eliminate the variable in the first place. 
	
		\begin{example}
		Here is a simple example which illustrates this technique. Consider the following LP program which we wish to standardise:

		\[
			\begin{array}{lrrrrrrr}
			\mathrm{Minimise}    &  x_1 & + & 3x_2 & + & 4x_3 &   &   \\
			\mathrm{subject\ to} &  x_1 & + & 2x_2 & + &  x_3 & = & 5 \\
	         			         & 2x_1 & + & 3x_2 & + &  x_3 & = & 6 \\
	         			         & \mathrlap{x_2, x_3\ge 0.} & & & & & &
			\end{array}
		\]
		
		Since $x_1$ is an unconstrained variable, use the first equation to substitute $x_1 = 5 - 2x_2 - x_3$ and convert the original LP problem into:
		
		\[
			\begin{array}{lrrrrr}
			\mathrm{Minimise}    &  x_2 & + & 3x_3 & + & 5 \\
			\mathrm{subject\ to} &  x_2 & + &  x_3 & = & 4 \\
	         			         & \mathrlap{x_2, x_3\ge 0,} & & & &
			\end{array}\\[0.2cm]
		\]
		
		which is now in standard form.
		
		\end{example}

\end{enumerate}

\subsection{The Fundamental Theorem of LP}

In this section we will develop the necessary background on basic solutions for linear programming problems that will allow us to state the Fundamental Theorem of Linear Programming (see Section 2.3 from \cite{ref7}).
\\

Let us therefore focus on the standard form given by \ref{standard-form-compact}. Recall that $\bm{A}\in\mathbb{R}^{m\times n}$, $\bm{b}\in\mathbb{R}_+^m$ and $\bm{c}\in\mathbb{R}^n$. In what follows, we will assume that the number of decision variables is always greater than or equal to the number of equations, i.e. $n\ge m$, as otherwise the system $\bm{Ax} = \bm{b}$ will be overdetermined. Also, we assume that a problem written in standard form has no redundant or inconsistent constraints, i.e. the rows of matrix $\bm{A}$ are linearly independent and hence $\text{rank}(\bm{A}) = m$.
\\

Now, consider the system of linear equations from the standard form \ref{standard-form-compact}:

\begin{equation} \label{Ax=b}
\bm{Ax} = \bm{b}\\[0.2cm]
\end{equation}

\noindent and suppose, for convenience, that we denote by $\bm{B}$ the $m\times m$ matrix formed by the first $m$ columns of matrix $\bm{A}$. Then it follows that the matrix $\bm{B}$ is non-singular, so the equation:

\[\bm{Bx_B} = \bm{b}\\[0.2cm]\]

\noindent is guaranteed to have a unique solution $\bm{x_B}$. By setting the first $m$ elements of the vector $\bm{x}$ to be equal to those of $\bm{x_B}$ and filling the remaining $n - m$ entries with $0$, i.e. $\bm{x} = [\bm{x_B}, \bm{0}]$, we have a solution for the linear system $\bm{Ax} = \bm{b}$. We are now in a position to formulate:

\begin{definition}
Given the set of $m$ simultaneous linear equations in $n$ unknowns \ref{Ax=b}, let $\bm{B}$ be any nonsingular $m\times m$ sub-matrix consisting of columns of $\bm{A}$. Then, if all $n - m$ entries of $\bm{x}$ not associated with columns of $\bm{B}$ are set equal to zero, the solution to the resulting set of equations is said to be a \textit{basic solution} to \ref{Ax=b} with respect to the basis $\bm{B}$. The entries of $\bm{x}$ associated with columns of $\bm{B}$ are called \textit{basic variables}.	
\end{definition}

We also need to note that, under the rank assumption for matrix $\bm{A}$, i.e. $\text{rank}(\bm{A}) = m$, the system of linear equations \ref{Ax=b} will always have at least one basic solution. However, basic variables in a basic solution can potentially be zero as well. We thus introduce the following:

\begin{definition}
	If one or more of the basic variables in a basic solution is equal to zero, that solution is called \textit{degenerate basic solution}.
\end{definition}

Until this point, we have not yet treated the non-negativity constraints of the decision variables, i.e. $\bm{x} \ge 0$. So let:

\begin{equation} \label{all-constraints}
\begin{array}{c}
\bm{Ax} = \bm{b} \\
\bm{x} \ge \bm{0}
\end{array}
\end{equation}

be the system of constraints for any LP problem in the standard form \ref{standard-form-compact}. We can now state similar definitions that apply when these restrictions are considered together, so let us present:

\begin{definition}
	A vector $\bm{x}$ satisfying \ref{all-constraints} is called \textit{feasible} for these constraints. A feasible solution to the constraints \ref{all-constraints} that is also basic is called \textit{basic feasible solution}; if this solution is also a degenerate basic solution, it is called a \textit{degenerate basic feasible solution}.
\end{definition} 

Before introducing the fundamental theorem of LP, we require one more definition to account for the optimality of solutions for an LP problem, so for this reason we present:

\begin{definition}
	Given an LP in standard form, a feasible solution to the constraints \ref{all-constraints} that achieves the optimal value of the objective function is called an \textit{optimal feasible solution}. If the solution is basic then it is an \textit{optimal basic feasible solution}.
\end{definition}

In this moment we have all the necessary concepts to formulate the main result from this section:

\begin{theorem} (Fundamental Theorem of Linear Programming) Given a linear program in standard form \ref{standard-form-compact}, where $\bm{A}$ is an $m\times n$ matrix with $\text{rank}(A) = m$,

\begin{enumerate}
	\item if there is a feasible solution, there is a basic feasible solution;
	\item if there is an optimal feasible solution, there is an optimal basic feasible solution.
\end{enumerate}
	
\end{theorem}

\subsection{Polyhedral Convex Sets}

In this section we will depart from the elementary properties of linear systems that we employed in the study of the fundamental theorem of linear programming and rather focus on a more intuitive geometrical interpretation. Concepts from the standard theory of convex sets will be used to give an alternative definition to the fundamental theorem, which will lead to a clearer geometric understanding of the result. As it will soon become apparent, there is an intimate link between the basic feasible solutions of the algebraic theory and extreme points of convex polygons in the geometric approach. We proceed with the following definitions and theorems, followed in the next section by a concrete example.

\begin{definition} Let $\bm{x} = (x_1, x_2, \ldots, x_n)$, $\bm{y} = (y_1, y_2, \ldots, y_n)\in\mathbb{R}^n$. Then the \textit{line segment} between $\bm{x}$ and $\bm{y}$ (including the endpoints) is given by:

\[\alpha\bm{x} + (1-\alpha)\bm{y},\ 0\le \alpha\le 1.\] 

\end{definition}

\begin{definition} Let $S\subset \mathbb{R}^n$. The set $S$ is said to be \textit{convex} if for all $\bm{x} = (x_1, x_2, \ldots, x_n)$, $\bm{y} = (y_1, y_2, \ldots, y_n)\in S$ it holds that:

\[\alpha\bm{x} + (1-\alpha)\bm{y},\ 0\le \alpha\le 1.\]

\end{definition}

\begin{definition} The set of points $(x_1, x_2, \ldots, x_n)\in\mathbb{R}^n$ satisfying an equation of the form

\[a_1x_1 + a_2x_2 + \ldots + a_nx_n = b\]

\noindent is said to be a \textit{hyperplane} of $\mathbb{R}^n$. The set of points $(x_1, x_2, \ldots, x_n)\in\mathbb{R}^n$ satisfying an inequality of the form

\[a_1x_1 + a_2x_2 + \ldots + a_nx_n \lessgtr b\]

\noindent is said to be a \textit{closed half-space} of $\mathbb{R}^n$.
\end{definition}

\begin{theorem} The constraint set of a canonical maximisation/minimisation linear programming problem is convex. Such a set is said to be a polyhedral convex set.	
\end{theorem} 

\begin{definition} Let $\bm{x} = (x_1, x_2, \ldots, x_n)\in\mathbb{R}^n$. The norm of $\bm{x}$, denoted $\Vert \bm{x} \Vert$, is given by:

\[\Vert \bm{x} \Vert = \sqrt{x_1^2 + x_2^2 + \ldots + x_n^2}.\]

\end{definition}

\begin{definition} Let $r\ge 0$. The set of points $\bm{x} = (x_1, x_2, \ldots, x_n)\in\mathbb{R}^n$ such that

\[\Vert \bm{x} \Vert \le r\]

\noindent is said to be the \textit{closed ball of radius} $r$ \textit{centred at the origin}.
\end{definition}

\begin{definition} A set $S\subset\mathbb{R}^n$ is said to be \textit{bounded} if there exists $r\ge 0$ such that every element of $S$ is contained in the closed ball or radius $r$ centred at the origin. A set $S\subset\mathbb{R}^n$ is said to be \textit{unbounded} if it is not bounded.
\end{definition}

\begin{definition} Let $S$ be a convex set in $\mathbb{R}^n$. A point $e\in S$ is said to be an extreme point of $S$ if there do not exist $\bm{x}, \bm{y}\in S$ and $\alpha\in (0, 1)$ such that

\[\textbf{e} = \alpha \textbf{x} + (1-\alpha) \textbf{y}.\]

\end{definition}

\begin{theorem}\label{T1} If the constraint set $S$ of a linear programming problem in standard form is bounded, then the maximum/minimum value of the objective function is attained at an extreme point of $S$.	
\end{theorem}

\begin{theorem}\label{T2} If the constraint set $S$ of a linear programming problem in standard form is unbounded, then there exists some $M\in\mathbb{R}$ such that the objective function $f$ satisfies $f(x_1, x_2, \ldots, x_n)\lessgtr M$ for all $(x_1, x_2, \ldots, x_n)\in S$, i.e. $f$ is bounded above/below (by $M$), them the maximum/minimum value of the objective function is attained at an extreme point of $S$.	
\end{theorem}

\subsection{A Geometric Approach}

Let us illustrate the above concepts with a typical example of resource allocation problem.

\begin{example} A furniture company manufactures sofas and armchairs. The production of one sofa requires 2 hours in the parts division and 1 hour on the assembly line of the company. The production of one armchair requires 1 hour in the parts division and 2 hours on the assembly line. The parts department of the company operates at most 8 hours per day, while the assembly division is active at most 10 hours per day. Knowing that the profit of selling one sofa is \pounds 30 and the profit of selling one armchair is \pounds 50, what are the optimal quantities of sofas and armchairs that the company should produce in order to maximise profits?	
\end{example}

We start by translating the problem in mathematical terms. Note that we are interested in the number of sofas and armchairs to be produced, so let us denote:

\[\begin{array}{lll}
x_1 & = & \text{\# of sofas per day;} \\ 
x_2 & = & \text{\# of armchairs per day.}
\end{array}
\]

According to the above definitions, $x_1$ and $x_2$ are the \textbf{decision variables} of our problem. Next, we wish to maximise the profits of the company, namely: 

\[f(x_1, x_2) = 30 x_1 + 50 x_2.\]

This is the \textbf{objective function} that we want to optimise. But we cannot have unlimited quantities and hence infinite profits, since the company is constrained by the availability of the parts and assembly divisions. We observe that 2 hours are spent in the parts division for one sofa and 1 hour is spent in the same division for one armchair. Together with the constraint of 8 hours per day in this department, we derive the constraint:

\[2x_1 + x_2 \le 8.\]

Similarly, we can derive the other constraint for the assembly line, which gives a second inequality:

\[x_1 + 2x_2 \le 10.\]

Finally, we also have the obvious constrains:

\[
\begin{array}{c}
x_1 \ge 0, \\x_2 \ge 0. 	
\end{array}
\]

since the company cannot produce negative quantities of sofas or armchairs. The above 4 constraints represent the \textbf{feasible set} that enforce the admissible production plans $x = (x_1, x_2)$. Note that, for example, $x = (20, 30)$ cannot belong to this feasible set, as there is not enough availability in both the parts and assembly divisions to satisfy that level of production. The final optimisation problem is a maximisation problem that can be formulated mathematically as follows:

\begin{figure}[H]
	\centerline{\includegraphics[scale=0.5]{pic1}}
\label{figure1}
\caption{Feasible set of points satisfying all constraints}
\label{figure1}
\end{figure}

\[
\begin{array}{ll}
\text{Maximise}   & f(x_1, x_2) = 30x_1 + 50x_2 \\
\text{subject to} & 2x_1 + x_2 \le 8 \\
	              & x_1 + 2x_2 \le 10 \\
	              & x_1 \ge 0 \\
	              & x_2 \ge 0.
\end{array}
\]

The set of points $(x_1, x_2)$ satisfying the above constraints is given by the shaded region from Figure \ref{figure1}.
\\

\begin{figure}[!htb]
\minipage{0.5\textwidth}
	\centerline{\includegraphics[scale=0.35]{pic2}}
\endminipage\hfill
\minipage{0.5\textwidth}
	\centerline{\includegraphics[scale=0.35]{pic3}}
\endminipage
\end{figure}

\vspace{0.2cm}

\begin{figure}[!htb]
\minipage{0.5\textwidth}
	\centerline{\includegraphics[scale=0.35]{pic4}}
\endminipage\hfill
\minipage{0.5\textwidth}
	\centerline{\includegraphics[scale=0.35]{pic5}}
\endminipage
\caption{All individual constraints}
\label{figure2}
\end{figure}

The region in Figure 1 above was obtained by intersecting the four 2-dimensional regions determined by the constraints $2x_1 + x_2 \le 8$, $x_1 + 2x_2 \le 10$, $x_1 \ge 0$, and $x_2 \ge 0$ (Figure \ref{figure2}). Therefore, the question boils down to determining which point in the shaded region in Figure 1 maximises $f(x_1, x_2)$.
\\

Since the constraint set is bounded, we can use Theorem \ref{T1} to deduce that the maximum value of $f(x, y)$ is attained at an extreme point in Figure 1, namely at either $(0, 0)$, $(4, 0)$, $(0,5)$, or $(2,4)$. A simple analysis shows that the maximum is attained at $(2, 4)$ where $f(2, 4) = 260$, and as a result the company should produce $2$ sofas and $4$ armchairs in order to maximise their profits (and get $\pounds 260$ in return).

\subsection{Complexity Limitations and Alternative Approaches}

Having seen the geometric method which allows us to reason about LP problems, we can already notice some of its limitations. In real-life LP problems, one would typically deal with potentially tens if not hundreds of variables. 
\\

For example, assuming $m$ constraints and $n$ decision variables, we will have at most

\[{m + n \choose n} = \frac {(m + n)!}{m!\cdot n!}\]

\noindent candidate extrema points to be tested. This clearly poses a problem and makes the visualisation of the constraint set impossible.
\\

Another more relevant problem concerns the cases when the constraint set is unbounded. In such scenarios, according to Theorem \ref{T2} we would also need to ensure that the objective function is bounded by above or below, depending on whether the LP problem at hand is a maximisation or a minimisation problem, respectively. This is largely problem dependant and can complicate the analysis significantly. One such example is given by the following LP problem:

\[
\begin{array}{ll}
\text{Maximise}   & f(x, y, z) = 2x + 3y + 4z \\
\text{subject to} &  y + 5x \le 10 \\
	              & 2y + 3z \le 15 \\
	              & x, y, z \ge 0
\end{array}
\]

Here there is no restriction on $x$ other than being positive, and as such for $x\to\infty$ we will have $f(x, y, z)\to\infty$, and hence the objective function is unbounded.
\\

Given these obvious limitations, several algorithms have been developed to address this rather exponential nature of linear programming problems. The famous \textit{simplex algorithm} is able to find optimal solutions to LP problems without testing a large number of candidate extrema points. It is also able to detect edge cases where the constraint sets are empty and the objective functions are unbounded (see \cite{ref6}). The basic idea is that the simplex method is capable of continually decreasing the value of the objective function by intelligently \textit{pivoting} from one feasible solution (i.e. one extreme point) to another. This reduces the search space considerably and leads to a more computationally-friendly algorithm. However, even with such optimisations in place, the simplex algorithm can still exhibit exponential time complexity, as shown by Klee-Minty (see Section 5.2 from \cite{ref6} for such an example).
\\

Cutting-plane algorithms such as the ellipsoid method have been found to have polynomial-time complexity in the size of the problem ($\mathcal{O}((m + n)^4)$), but it was discovered that in practice they do not perform any better than the simplex algorithm. 
\\

Currently, the state-of-the-art algorithm that is being used by major LP software packages relies on  the interior-point method belonging to Karmarkar which has total complexity of the order $\mathcal{O}(nm^2\log(n/\varepsilon))$, where $\varepsilon$ is a tolerance parameter (see also \cite{ref6}). This method was found to behave much better in practice than the simplex algorithm and it is also reasonably effective when applied to large-scale LP problems.

\subsection{Duality Theory}

In this last subsection on linear programming, we briefly discuss the duality relationship that holds for certain pairs of LP problems. To this end, we will deviate from the standard form given by \ref{standard-form}, as duality arises naturally from the symmetry of LP problems expressed only in terms of inequalities. We now define duality through the following pair of programs:

\begin{equation}
\label{duals}
\begin{aligned}
\begin{array}{ll}
\text{Maximise}   & \bm{c}^\top \bm{x} \\
\text{subject to} & \bm{Ax} \le \bm{b} \\
	              & \bm{x} \ge \bm{0}
\end{array} \qquad & \qquad
\begin{array}{ll}
\text{Minimise}   & \bm{b}^\top \bm{y} \\
\text{subject to} & \bm{A}^\top\bm{y} \ge \bm{c} \\
	              & \bm{y} \ge \bm{0}
\end{array} 
\end{aligned},\\[0.2cm]
\end{equation}

\noindent where $\bm{A}$ is an $m\times n$ matrix, $\bm{x}$ is an $n$-dimensional column vector and $\bm{c}, \bm{y}$ are $m$-dimensional column vectors. The symmetry of the above LP problems allows us to define a \textit{duality} relationship between the two programs: one of them will be called \textit{primal}, while its counterpart will be named \textit{dual}. It is also important to note that the roles of the primal and dual problems can be swapped. We are now able to state the fundamental theorems on duality:


\begin{theorem} (Weak Duality). If $\bm{x}$ and $\bm{y}$ are feasible solution for the LP problems defined by \ref{duals}, then the following inequality holds: $\bm{c}^\top \bm{x} \le \bm{b}^\top\bm{y}$. 
	
\end{theorem}

\begin{proof}
	Since $\bm{x}\ge 0$ we have: $\bm{c}^\top \bm{x}\le \left(\bm{A}^\top \bm{y}\right)^\top \bm{x} = \bm{y}^\top \bm{A}\bm{x} \le \bm{y}^\top \bm{b} = \bm{b}^\top \bm{y}.$
\end{proof}

\vspace{0.2cm}

\begin{theorem} (Strong Duality). If either of the problems defined in \ref{duals} has a finite optimal solution, so does the other, and the corresponding values of the objective functions coincide. If either problem has an unbounded objective, the other problem has no feasible solution as well.
\end{theorem}

\begin{proof}
	The interested reader can refer to Section 4.2 from \cite{ref7}.	
\end{proof}

\vspace{0.2cm}

We conclude this brief tour of linear programming with some remarks that are connected to the duality relationship presented above. In the view of Theorem \ref{min-max} presented in an earlier section, one can simplify the computational framework for LP problems by treating all such problems as being either minimisation problems or maximisations problems. Usually, most optimisations problems are formulated as minimisations problems and a convention is followed as to whether the column vector that defines the inequality constraints represents the lower bounds or the upper bounds of those constraints. Generally, the preferred way is to define the constraints as: $\bm{G}\bm{x} \le \bm{h}$, so that $\bm{h}$ represents the column vector of right hand side terms. 

\section{CVXOPT Framework} \label{cvxopt-section}

In order to implement the linear programming algorithms that decide consistency of function and derivative information (which will be described in a later section), we had to choose a library that is specifically tailored for dealing with LP problems. Previous experience with GNU's GLPK package \cite{ref9}, which uses a rather obfuscated syntax with many intricate constructs, made us look for an alternative package in the open-source landscape. Given that Python is a well-established programming language, which is known for its ease of use which allows fast and convenient code development, we decided to narrow down our search to linear programming solvers available in Python.
\\

After a quick search for such libraries, we discovered that the CVXOPT library \cite{ref8} will be suitable for the purposes of this project, since it combines both decent performance and convenient formulation of the LP problems. In addition to handling problems with linear objectives, CVXOPT also features support for more general convex optimisation problems which have non-linear objective functions [citation needed].

\subsection{Formulation of LP Problems in CVXOPT}

The API for specifying the constraints and the objective function of a linear programming problem is very simple and has the following function signature (which is simplified here for presentation purposes) [citation needed]:
\[\verb|cvxopt.solvers.lp(c, G, h[, A, b[, solver]])|,\]

\noindent and is designed for solving the following minimisation problem:

\begin{equation}
\label{LP-cvxopt}
\begin{array}{ll}
\text{Minimise}   & \bm{c}^\top \bm{x} \\
\text{subject to} & \bm{Gx} \le \bm{h} \\
				  & \bm{Ax} = \bm{b} \\
\end{array}.\\[0.2cm]
\end{equation}

\noindent Note that the parameters \lit*{c}, \lit*{G}, \lit*{h} are compulsory, while the remaining ones are optional. Depending on whether the LP problem at hand features equality constraints or not, one would need to specify both the following parameters \lit*{A} and \lit*{b}. In addition, external solvers such as \lit{glpk} or \lit{mosek} (if installed) can be used instead of the default one used internally, by means of the argument \lit*{solver}. For completeness, assuming that there are $n$ decision variables, $p$ inequality constraints and $m - p$ inequality constraints, it holds that: $\bm{c}\in\mathbb{R}^n$, $\bm{G}\in\mathbb{R}^{p\times n}$, $\bm{h}\in\mathbb{R}^p$, $\bm{A}\in\mathbb{R}^{(m - p)\times n}$ and $\bm{b}\in\mathbb{R}^{m - p}$.

\subsection{Simple Example}

We can now illustrate the strengths of the CVXOPT framework with a simple example. Given the following optimisation problem:

\[
\begin{array}{ll}
\text{Minimise}   & f(x_1, x_2) = 2x_1 + x_2 \\
\text{subject to} & -x_1 + x_2 \le 1 \\
	              & x_1 + x_2 \ge 2 \\
	              & x_2 \ge 0 \\
	              & x_1 - 2x_2 \le 4,
\end{array}\\[0.2cm]
\]

\noindent we ask for the optimal solution vector $\bm{x} = (x_1, x_2)\in\mathbb{R}^2$ which minimises the given linear objective function. Since we are only dealing with inequalities constraints, we simply need to construct the vector $\bm{c}\in\mathbb{R}^2$, the coefficient matrix $\bm{G}\in\mathbb{R}^{2\times 4}$ and the column vector of right hand sides $\bm{h}\in\mathbb{R}^4$. By transforming the above optimisation problem into the form \ref{LP-cvxopt}, we have:
\[
\begin{aligned}
\bm{c} = \begin{bmatrix*}[r]
 2 \\ 1	
\end{bmatrix*}
 \quad && \quad
\bm{G} = \begin{bmatrix*}[r]
 -1 & 1 \\ -1 & -1 \\ 0 & -1 \\ 1 & 2	
\end{bmatrix*}
 \quad && \quad
\bm{h} = \begin{bmatrix*}[r]
 1 \\ -2 \\ 0 \\ 4
\end{bmatrix*}
\end{aligned},
\]

\noindent which we can now input in the interactive Python interpreter, followed by the call to \lit*{solvers.lp} to retrieve the optimal solution to the given problem:

\begin{lstlisting}
>>> from cvxopt import matrix, solvers
>>> c = matrix([2.0, 1.0])
>>> G = matrix([[-1.0, -1.0, 0.0, 1.0], [1.0, -1.0, -1.0, -2.0]])
>>> h = matrix([1.0, -2.0, 0.0, 4.0])
>>> sol = solvers.lp(c, G, h)

     pcost       dcost       gap    pres   dres   k/t
 0:  2.6471e+00 -7.0588e-01  2e+01  8e-01  2e+00  1e+00
 1:  3.0726e+00  2.8437e+00  1e+00  1e-01  2e-01  3e-01
 2:  2.4891e+00  2.4808e+00  1e-01  1e-02  2e-02  5e-02
 3:  2.4999e+00  2.4998e+00  1e-03  1e-04  2e-04  5e-04
 4:  2.5000e+00  2.5000e+00  1e-05  1e-06  2e-06  5e-06
 5:  2.5000e+00  2.5000e+00  1e-07  1e-08  2e-08  5e-08

>>> print(sol[`x'])
[ 5.00e-01]
[ 1.50e+00]

\end{lstlisting}

\section{Domain Theory}

We introduce basic domain theory knowledge that will become useful when dealing with the problem of consistency in a later section. We use \cite{ref5} and \cite{ref10} when referencing these results, unless otherwise stated.

\subsection{Introduction}

Domain theory was introduced in 1970 by Dana Scott as a mathematical theory of programming languages.  The basic idea of domain theory is to provide better and better approximations to an object by means of simple recursion. It has applications in the field of computer science, e.g. solving canonically fixed point equations or recursive equations of procedures and data structures.

\subsection{Main Definitions and Examples}

\begin{definition} A \textit{partial order} (or a \textit{partially ordered set} or \textit{poset}) $(D, \le)$ is a set $D$ together with a binary relation $\le$ which is:

\begin{enumerate}
	\item reflexive: $a\le a$,
	\item anti-symmetric: $a\le b \wedge b\le a \implies a = b$, and
	\item transitive: $a\le b \wedge b\le c \implies a\le c$.
\end{enumerate}

Most often, the binary relation of a partial order is written as $\sqsubseteq$. Then $a\sqsubseteq b$ can be interpreted as $a$ having less information than $b$.
\end{definition} 

\begin{definition} A subset $A$ of an ordered set $(P, \sqsubseteq)$ is an \textit{upper set} if $x\in A\implies y\in A$, for all $y\sqsupseteq x$. We denote by $\uparrow A$ the set of all elements above some element of $A$. For convenience, we abbreviate $\uparrow \{x\}$ as $\uparrow x$. The dual notions are lower set and $\downarrow A$ \cite{ref11}.
\end{definition} 

\begin{definition} A non-empty subset $A\subset P$ is said to be \textit{directed} if for any pair of elements $a, b\in A$ there exists $c\in A$ such that $a\sqsubseteq c$ and $b\sqsubseteq c$.
\end{definition} 

\begin{definition} A \textit{directed complete partial order} (\textit{dcpo}) or a domain is a partial order in which every directed subset has a least upper bound (lub). A dcpo is said to be \textit{pointed} if it has a least element which is denoted by $\bot$ and is called bottom.
\end{definition}

\begin{definition} Given two elements $a$ and $b$ of a dcpo we say that $a$ is \textit{way-below} or \textit{approximates} $b$, denoted by $a\ll b$, if for every directed subset $A$ with $b\sqsubseteq \bigsqcup A$ there exists $c\in A$ with $a\sqsubseteq c$.
\end{definition} 

\begin{definition} A basis of a domain $D$ is a subset $B\subset D$ such that for every element $x\in D$ of the domain the set $B_x = \{y\in B\ |\ y \ll x\}$ of elements in the basis way-below $x$ is directed with $x = \bigsqcup B_x$. Also, a dcpo with a (countable) basis is said to be an \textit{($\omega$)-continuous} domain. 

\end{definition} 

\begin{definition} A function $f : D \to E$ between dcpo's is said to be \textit{Scott-continuous} if and only if it is \textit{monotone} (i.e. $a\sqsubseteq b \implies f(a) \sqsubseteq f(b)$) and preserves' lub's of directed sets i.e. for any directed $A\subseteq D$, we have $f\left(\bigsqcup_{a\in A} a\right) = \bigsqcup_{a\in A} f(a)$. Moreover, if $D$ is an $\omega$-continuous dcpo, then $f$ is continuous if and only if it is monotone and preserves the lub's of increasing sequences (i.e. $f\left(\bigsqcup_{i\in\omega} x_i\right) = \bigsqcup_{i\in\omega} f(x_i)$, for any increasing $(x_i)_{i\in\omega}$).
\end{definition} 

\begin{definition} Let $D$ be a dcpo. A subset $A$ is called \textit{(Scott-)closed} if it is a lower set and is closed under suprema of directed subsets. Complements of closed sets are called \textit{(Scott-)open}; they are the elements of $\omega_D$, the \textit{Scott-topology} on $D$ \cite{ref11}.
\end{definition} 

\begin{example} \cite{ref5}
	The interval domain $\mathbf{I}[0, 1]^n$ of the unit box $[0, 1]^n \subseteq \mathbb{R}^n$ is the set of all non-empty $n$-dimensional sub-rectangles in $[0, 1]^n$ ordered by reverse inclusion. A basic Scott open set is given, for every open subset $O$ of $\mathbb{R}^n$, by the collection of all rectangles contained in $O$. 
	\\
	
	The map $x \mapsto \{x\} : [0, 1]^n \to \mathbf{I}[0, 1]^n$ is an embedding onto the set of maximal elements of $\mathbf{I}[0, 1]^n$ . Every maximal element $\{x\}$ can be obtained as the least upper bound (lub) of an increasing chain of elements, that is a shrinking, nested sequence of sub-rectangles, each containing $\{x\}$ in its interior and thereby giving an approximation to $\{x\}$ or equivalently to $x$. The set of sub-rectangles with rational coordinates provides a countable basis. One can similarly define, for example, the interval domain $\mathbf{I}\mathbb{R}^n$ of $\mathbb{R}^n$.
\end{example}






