I have an algorithm, but no closed formula for you.
Algorithm
Compute the convex hull of your points. The ellipsis you search for has to touch the corners of that hull in several points if it is to be the smallest one. Now iterate over all 4-element subsets of those corners, and solve the following equation:
$\begin{pmatrix} x_1^2 & y_1^2 & x_1 & y_1 & 1 \\ x_2^2 & y_2^2 & x_2 & y_2 & 1 \\ x_3^2 & y_3^2 & x_3 & y_3 & 1 \\ x_4^2 & y_4^2 & x_4 & y_4 & 1 \end{pmatrix}\cdot\begin{pmatrix} a\\b\\c\\d\\e \end{pmatrix}=0 $
The result won't be unique, but instead a one-dimensional space of solutions unless you hit a degenerate situation. Take any non-zero vector from that space to describe your conic secion. You may e.g. add a fifth equation $a=1$ as long as you ensure that you will recover from any degeneracies (i.e. unsolvable systems due to linear dependency).
The conic section described by the parameters you found will already be “parallel” to the coordinate axes, and it can be written as
$ ax^2 + by^2 + cx + dy + e = 0 $
Check whether this conic section is actually an ellipsis, and if it contains all other points of the convex hull. If it is and does, compute its area, and compare that against the best solution found so far.
Details
To perform all the computations you need, let's turn the above formula into something easier to manage.
\begin{align*} 0 &= ax^2 + by^2 + cx + dy + e \\ &= a\left(x+\frac{c}{2a}\right)^2 + b\left(y+\frac{d}{2b}\right)^2 - \left(\frac{c^2}{4a} + \frac{d^2}{4b} - e\right) \\ m_x &= -\frac{c}{2a} \\ m_y &= -\frac{d}{2b} \\ s &= \left(\frac{c^2}{4a} + \frac{d^2}{4b} - e\right) \end{align*}
$(m_x, m_y)$ describes the center of symmetry for your conic.
If it is to be an ellipsis, then $a$ and $b$ have to have the same non-zero sign (i.e. $ab>0$), in which case $s$ will have the same sign as well. Now divide your whole equation by $s$ to obtain
\begin{align*} \frac as(x-m_x)^2 + \frac bs(y-m_y)^2 + \frac ss &= 0 \\ \frac as(x-m_x)^2 + \frac bs(y-m_y)^2 &= 1 \\ a'(x-m_x)^2 + b'(y-m_y)^2 &= 1 \end{align*}
with
$a' = \frac as \quad b' = \frac bs$
denoting the inverse of the square lengths of the radii of your ellipsis.
You can turn the equality into an inequality $\ldots < 1$ to describe the interior of the ellipsis, and check whether the other points of the convex hull fall inside. You can use the term $a'\cdot b'$ to compare sizes, as it is indirectly proportional to the squared area of your ellipsis.
If you need the actual radii, you can compute them as
\begin{align*} r_x &= \frac1{\sqrt{a'}} = \sqrt{\frac sa} \\ r_y &= \frac1{\sqrt{b'}} = \sqrt{\frac sb} \\ 1 &= \left(\frac{x-m_x}{r_x}\right)^2 + \left(\frac{y-m_y}{r_y}\right)^2 \end{align*}
Alternatives
For users who read this question while investigating a different task:
- For arbitrary center and fixed orientation, we used four points to define the conic, as seen above.
- For arbitrary center and arbitrary orientation, we need five points, as five points uniquely define a general conic. Add a $xy$ monomial to the equation.
- For a fixed center like the origin, and fixed orientation, use two points, as $c=d=0$.
- For fixed center and arbitrary orientation, use three points.