Roundest ellipse with specified tangents
I found a way to do this construction in a purely geometrical way, without coordinates. Let $A$, $B$ be the given tangency points (called $\mathbf p_1$ and $\mathbf p_2$ in the question) and $AV$, $BV$ the corresponding tangents (the case when tangents are parallel can be dealt with in a simpler way, see below).
If $M$ is the midpoint of $AB$, for a well-known property of the ellipse line $VM$ is a diameter, and we are free to choose on it a point $D$, to construct the ellipse touching the given tangents at $A$, $B$ and passing through $D$. The construction is quite simple, because the line through $D$ parallel to $AB$, meeting line $VA$ at $E$, is another tangent: if $F$ is the midpoint of $AD$ it follows that $EF$ is also a diameter, and the intersection $C$ of $VM$ and $EF$ is thus the center of the ellipse (see diagram). With the center and three points, the ellipse can be easily constructed.
But of course we must choose point $D$ such that the eccentricity of the ellipse be minimum. To that end, draw line $GC$ parallel to $AB$ and notice that lines $CD$ and $CG$ form a pair of conjugate diameters. Apollonius' equations connect semidiameters $CD$ and $CG$ with standard semiaxes $a$ and $b$ of the ellipse: $$ CD^2+CG^2=a^2+b^2,\quad CD\cdot CG\sin\theta = ab, $$ where $\theta=\angle VMA$. From there it is not difficult to find an expression for the eccentricity $\epsilon$ as a function of $t=CD/CG$: $$ \epsilon^2={2\over1+1/\xi}, \quad\hbox{where}\quad \xi=\sqrt{1-\left({2t\over1+t^2}\right)^2\sin^2\theta}. $$ It follows that
the minimum value of the eccentricity, $\epsilon=\sqrt{2\over1+1/|\cos\theta|}$, is obtained when $CD=CG$.
To exploit that condition, notice that $(AM/CG)^2+(CM/CD)^2=1$ and if $CD=CG$ this is equivalent to $$ AM^2+CM^2=CD^2. $$ Playing around with similar triangles one also finds $CD/MC=1+VM/MD$: combining those equalities we finally get $$ MD={AM^2+AM\sqrt{AM^2+VM^2}\over VM}. $$ From that, one can construct point $D$ and then find center $C$ as explained above.
If tangent lines are parallel, $A$ and $B$ are the endpoints of a diameter and center $C$ is their midpoint. Just construct then the conjugate diameter through $C$ parallel to the tangents and take on it point $D$ such that $CD=AC$.
At the moment I can offer an analytical solution, which is quite cumbersome but manageable. I hope someone will be able to give a simpler solution.
To reduce the number of parameters, we can set up a coordinate system such that $\mathbf{p}_1=(\alpha,0)$, $\mathbf{p}_2=(0,\beta)$, $\mathbf{v}_1=(0,1)$, and define then $m=v_{2y}/v_{2x}$.
Let's start with the generic equation for a conic section, where we can set to $1$ the coefficient of $x^2$ because we are dealing with an ellipse: $$ x^2+By^2+Cxy+Dx+Ey+F=0. $$
We have four conditions: the conic passes through $\mathbf{p}_1$, $\mathbf{p}_2$ and is there tangent to $\mathbf{v}_1$, $\mathbf{v}_2$. These conditions translate into four equations: $$ \begin{align} \cases{ \alpha^2+\alpha D+F=0\\ \beta^2B+\beta E+F=0\\ 2m\beta B+\beta C+D+mE=0\\ \alpha C+E=0 } \end{align} $$ From there, one can find expressions for $B$, $C$, $D$, $E$ as a function of $F$. Plugging these into the formulas for the semiaxes $a$ and $b$ we can then compute the eccentricity $\epsilon=\sqrt{1-b^2/a^2}$. The smaller the eccentricity, the rounder the ellipse, so we must find the minimum of $\epsilon$ as a function of $F$.
I computed $d\epsilon/dF$ with Mathematica and found that it vanishes for $$ F=-\frac{\alpha ^2 \beta \left(\alpha ^2 \beta +\beta ^3-2 \alpha ^2 \beta m^2+2 \alpha ^3 m\right)}{\beta ^2 \left(\alpha ^2+\beta ^2\right)+2 m^2 \left(\alpha ^4+2 \alpha ^2 \beta ^2\right)+2 \alpha \beta m \left(\alpha ^2+2 \beta ^2\right)}. $$
Once $F$ is known, one can compute the other coefficients and find the equation of the "roundest" ellipse: $$ \begin{align} \cases{ B=\frac{\displaystyle\alpha ^2 \left(\beta ^2+\alpha ^2 \left(2 m^2+1\right)+2 \alpha \beta m\right)}{\displaystyle\beta ^2\left(\alpha ^2+\beta ^2\right)+2 m^2 \left(\alpha ^4+2 \alpha ^2 \beta ^2\right)+2 \alpha \beta m \left(\alpha ^2+2 \beta ^2\right)}\\ \\ C=\frac{\displaystyle2 \alpha ^2 m \left(-\alpha ^2+\beta ^2+2 \alpha \beta m\right)}{\displaystyle\beta ^2 \left(\alpha ^2+\beta ^2\right)+2 m^2 \left(\alpha ^4+2 \alpha ^2 \beta ^2\right)+2 \alpha \beta m \left(\alpha ^2+2 \beta ^2\right)}\\ \\ D=-\frac{\displaystyle2 \alpha ^2 m \left(2 \beta ^3+m \left(\alpha ^3+3 \alpha \beta ^2\right)\right)}{\displaystyle\beta ^2 \left(\alpha ^2+\beta ^2\right)+2 m^2 \left(\alpha ^4+2 \alpha ^2 \beta ^2\right)+2 \alpha \beta m \left(\alpha ^2+2 \beta ^2\right)}\\ \\ E=-\frac{\displaystyle2 \alpha ^3 m \left(-\alpha ^2+\beta ^2+2 \alpha \beta m\right)}{\displaystyle\beta ^2 \left(\alpha ^2+\beta ^2\right)+2 m^2 \left(\alpha ^4+2 \alpha ^2 \beta ^2\right)+2 \alpha \beta m \left(\alpha ^2+2 \beta ^2\right)} } \end{align} $$
A check with GeoGebra confirms that this is indeed the value of $F$ giving a minimum eccentricity.
EDIT.
Following the suggestion by Rahul, we can get a much simpler result. Choose a coordinate system such that $\mathbf{p}_1=(\alpha,0)$, $\mathbf{p}_2=(-\alpha,0)$, and define $m=v_{1y}/v_{1x}$, $n=v_{2y}/v_{2x}$.
With those choices, the equation of the ellipse can be written as: $$ x^2+By^2-{m+n\over mn}xy+\alpha{m-n\over mn}y-\alpha^2=0. $$ The eccentricity $\epsilon$ is then a function of $B$, and $d\epsilon/dB=0$ for $$ B=1+\frac{(m+n)^2}{2 m^2 n^2}. $$ In summary, the equation of the roundest ellipse turns out to be: $$ x^2+\left(1+\frac{(m+n)^2}{2 m^2 n^2}\right)y^2-{m+n\over mn}xy+\alpha{m-n\over mn}y-\alpha^2=0. $$
Added by Rahul:
Things become even simpler if we use the slope of the normals, $\mu=-1/m$ and $\nu=-1/n$: $$ x^2+\left(1+\tfrac12(\mu+\nu)^2\right)y^2+(\mu+\nu)xy+\alpha(\mu-\nu)y-\alpha^2=0. $$
In an answer to a previous not-quite-duplicate question, achille hui has given an elegant coordinate-free solution, which I summarize below.
Choose the origin at the intersection of the two desired tangent lines, so that $\mathbf p_1$ and $\mathbf p_2$ are interpreted as vectors from the intersection to the two input points. Let $\mathbf q_1,\mathbf q_2$ be the dual basis for $\mathbf p_1,\mathbf p_2$, that is, $\mathbf p_1\cdot \mathbf q_1 = \mathbf p_2\cdot \mathbf q_2 = 1$ and $\mathbf p_1\cdot \mathbf q_2 = \mathbf p_2\cdot \mathbf q_1 = 0$. The equation for the minimum-eccentricity ellipse is $$(\mathbf{q}_1\cdot \mathbf{x} - 1)^2 + (\mathbf{q}_2\cdot \mathbf{x} - 1)^2 + 2\alpha(\mathbf{q}_1\cdot \mathbf{x})(\mathbf{q}_2\cdot \mathbf{x}) = 1,$$ where $$\alpha = \frac{2\mathbf{p}_1\cdot\mathbf{p}_2}{\|\mathbf{p}_1\|^2 + \|\mathbf{p}_2\|^2}.$$
For the derivation, see the original post.
Least eccentric ellipses for geometric Hermite interpolation
John C. Femiani, Chia-Yuan Chuang, Anshuman Razdan
Computer Aided Geometric Design 29 (2012) 141–149