Hint.
Given an ellipsoid and a line
$$
x^T A x = 1\\
x=x_0+\lambda v
$$
relaxing the $x_0$ definition as $\bar x$ we have
$$
(\bar x + \lambda v)^TA(\bar x+\lambda v) = 1\Rightarrow\bar x^T A \bar x+2\lambda \bar x A v +\lambda^2v^T A v = 1
$$
solving for $\lambda$
$$
\lambda = \frac{-2\bar x A v\pm\sqrt{4(\bar x A v)^2-4(v^T A v)(\bar x^T A \bar x-1)}}{v^T A v}
$$
at tangency we have
$$
(\bar x A v)^2-(v^T A v)(\bar x^T A \bar x-1)=0
$$
now choosing $\bar x$ such that $\bar x^T A \bar x = 1$ we have the condition
$$
\bar x A v = 0
$$
so $\bar x$ and $v$ should be orthogonal regarding $A$. Now the distances are computed as distances between $\bar x$ and the line $x=x_0+\lambda v$
NOTE
The intersection of the line $x = x_0+\lambda v$ with the plane $ x A v = 0$ is at $\lambda^* = -\frac{x_0^T A v}{v^T A v}$ or $x^* = x_0-\frac{x_0^T A v}{v^T A v}v$
This now is a planar problem that can be worked out in a plane as $O x y$. To do that we can rotate the plane $x^T A v$ around a convenient axis as for instance the axis formed by the line $\frac{v_x}{a^2}x + \frac{v_y}{b^2}y=0$ with direction $r = (-\frac{v_y}{b^2},\frac{v_x}{a^2},0)$ by an angle given by $\theta = \arccos\left(\frac{v_x}{||v||}\right)$. This can be done using the Rodrigues rotation formula (see). Now in the $Oxy$ plane we can compute the distance between the projected $x^*$ point and the projected ellipse.
NOTE
Included a MATHEMATICA script showing the adopted solution main lines
parms = {vx -> 2, vy -> 2, vz -> 1, a -> 8, b -> 4, c -> 1, x0 -> 8, y0 -> 0, z0 -> 2};
r = 11;
v = {vx, vy, vz}/Norm[{vx, vy, vz}];
zaxis = {0, 0, 1};
rpt = 0.1;
A = DiagonalMatrix[{1/a^2, 1/b^2, 1/c^2}];
p0 = {x0, y0, z0};
p = {x, y, z};
plane = p.A.v;
ellipsoid = p.A.p - 1;
past = p0 - (p0.A.v)/(v.A.v) v;
line = p0 + lambda v;
line0 = line /. parms;
linerot = plane /. {z -> 0};
axisrot = {-D[linerot, y], D[linerot, x], 0};
vrot = axisrot/Norm[axisrot];
theta = -ArcCos[zaxis.normal/Norm[normal]] /. parms;
R = RotationMatrix[theta, vrot] /. parms;
normal = Grad[plane, p];
rotellipsoid = p.R.A.Transpose[R].p - 1 /. parms // FullSimplify;
ellipsexy = rotellipsoid /. {z -> 0};
rotline = R.line /. parms;
vlrot = R.v;
p0rot = R.p0;
p0rot0 = p0rot /. parms;
pastrot = R.past;
pastrot0 = pastrot /. parms;
and the graphics.
h = plane /. parms;
g = ellipsoid /. parms;
gr3 = ParametricPlot3D[line0, {lambda, -r, r}];
gr23 = ContourPlot3D[{h == 0, g == 0}, {x, -r, r}, {y, -r, r}, {z, -r, r},
MeshFunctions -> {Function[{x, y, z, f}, h - g]},
MeshStyle -> {{Thick, Blue}}, Mesh -> {{0}},
ContourStyle -> Directive[Orange, Opacity[0.2], Specularity[White, 30]],
PlotPoints -> 60, SphericalRegion -> True];
grp0 = Graphics3D[{Red, Sphere[(p0 /. parms), rpt], PointSize[0.08]}];
grpast = Graphics3D[{Black, Sphere[(past /. parms), rpt], PointSize[0.08]}];
gr4 = ParametricPlot3D[lambda (vrot /. parms), {lambda, -r, r}];
h = z;
g = rotellipsoid;
gr34 = ContourPlot3D[{h == 0, g == 0}, {x, -r, r}, {y, -r, r}, {z, -r, r},
MeshFunctions -> {Function[{x, y, z, f}, h - g]},
MeshStyle -> {{Thick, Blue}}, Mesh -> {{0}},
ContourStyle -> Directive[Pink, Opacity[0.2], Specularity[White, 30]],
PlotPoints -> 60, SphericalRegion -> True];
gr3rot = ParametricPlot3D[rotline, {lambda, -r, r}];
grp0rot = Graphics3D[{Red, Sphere[(p0rot0 /. parms), rpt], PointSize[0.08]}];
grpastrot = Graphics3D[{Black, Sphere[(pastrot0 /. parms), rpt], PointSize[0.08]}];
Show[grp0, gr23, gr3, grpast, PlotRange -> {{-r, r}, {-r, r}, {-r, r}}]
Show[gr34, gr3rot, gr4, gr3rot, grp0rot, grpastrot, PlotRange -> {{-r, r}, {-r, r}, {-r, r}}]
and the distance calculations
rotline = p0rot + lambda vlrot;
L = (p - rotline).(p - rotline) + mu ellipsexy /. parms;
vars = Join[p, {lambda, mu}];
grad = Grad[L, vars];
sols = Solve[grad == 0, vars] // N;
pts = Table[(rotline /. parms /. sols[[k]]), {k, 1, Length[sols]}];
ptsz0 = Table[({x, y, 0} /. sols[[k]]), {k, 1, Length[sols]}];
grpts = Table[Graphics3D[{Blue, Sphere[pts[[k]], rpt]}], {k, 1, Length[sols]}];
grptsz0 = Table[Graphics3D[{Blue, Sphere[ptsz0[[k]], rpt]}], {k, 1, Length[sols]}];
grdists = Table[ParametricPlot3D[
pts[[k]] mu + ptsz0[[k]] (1 - mu), {mu, 0, 1},
PlotStyle -> Dashed], {k, 1, 4}];
Show[gr34, gr3rot, gr4, gr3rot,grpts, grptsz0, grdists, PlotRange -> {{-r, r}, {-r, r}, {-r, r}}]
Best Answer
[Collecting the discussion from the comments into an answer.]
You can certainly expand the expression that you’ve got and collect terms to get a quadratic equation in $\lambda$: Setting $\mathbf p = \mathbf o-\mathbf c$ and $\mathbf Q=\mathbf R^T\mathbf A\mathbf R$, you have $$(\lambda\mathbf l+\mathbf p)^T\mathbf Q(\lambda\mathbf l+\mathbf p) = (\mathbf l^T\mathbf Q\mathbf l)\lambda^2+2(\mathbf l^T\mathbf Q\mathbf p)\lambda +\mathbf p^T\mathbf Q\mathbf p=1,$$ which has real solutions when $(\mathbf l^T\mathbf Q\mathbf p)^2\ge(\mathbf l^T\mathbf Q\mathbf l)(\mathbf p^T\mathbf Q\mathbf p-1)$.
On the other hand, you’re essentially constructing the ellipsoid by scaling, rotating and translating the unit sphere. Since you’re not interested in the actual points of intersection, you could apply the inverse transformation to the line and then use the point-line distance formula to determine the line’s distance from the origin. If this is no greater than one, then there is at least one intersection.
If $\mathbf x_1$ and $\mathbf x_2$ are two points on a line, then the distance of the point $\mathbf x_0$ from this line is given by the formula (from Wolfram MathWorld) $${\lvert (\mathbf x_1-\mathbf x_0)\times(\mathbf x_2-\mathbf x_0)\rvert \over \lvert\mathbf x_2-\mathbf x_1\rvert}.$$ If $\mathbf x_0$ is the origin, this reduces to $${\lvert \mathbf x_2\times\mathbf x_1\rvert \over \lvert \mathbf x_2-\mathbf x_1\rvert}.\tag1$$ Now, since in your question you have $\mathbf R^T\mathbf A\mathbf R$ for the matrix of the ellipse, the matrix $\mathbf R$ already represents the inverse rotation. That is, it aligns the ellipsoid with the coordinate axes. (If instead you meant for it to represent a rotation that maps the standard basis onto the principal axes, then the ellipsoid’s matrix should be $\mathbf R\mathbf A\mathbf R^T$ instead and the inverse map below should use $\mathbf R^{-1}=\mathbf R^T$ instead of $\mathbf R$.) Setting $\mathbf S=\mathbf A^{1/2} = \operatorname{diag}(1/a,1/b,1/c)$, the required inverse map is therefore $\mathbf x\mapsto \mathbf S\mathbf R(\mathbf x-\mathbf c)$. A convenient choice for the two points on the line are the images of $\mathbf o$ and $\mathbf o+\mathbf l$. Plugging them into (1), we get $${\lvert \mathbf S\mathbf R(\mathbf o-\mathbf c)\times \mathbf S\mathbf R(\mathbf o+\mathbf l-\mathbf c) \rvert \over \lvert \mathbf S\mathbf R(\mathbf o-\mathbf c)-\mathbf S\mathbf R(\mathbf o+\mathbf l-\mathbf c) \rvert} = {\lvert \mathbf S\mathbf R(\mathbf o-\mathbf c)\times\mathbf S\mathbf R\,\mathbf l \rvert \over \lvert \mathbf S\mathbf R\,\mathbf l\rvert}.$$ Since we’re comparing a ratio of nonnegative numbers to $1$, we can square and rearrange to get $$\lvert\mathbf S\mathbf R(\mathbf o-\mathbf c)\times\mathbf S\mathbf R\,\mathbf l \rvert^2 \le \lvert \mathbf S\mathbf R\,\mathbf l\rvert^2.$$ This looks like it could be more efficient to compute that the quadratic equation’s discriminant.
We can even go a bit farther. Using $\lvert\mathbf u\times\mathbf v\rvert^2=\lvert\mathbf u\rvert^2\lvert\mathbf v\rvert^2-(\mathbf u\cdot\mathbf v)^2$, $$\lvert\mathbf S\mathbf R(\mathbf o-\mathbf c)\times\mathbf S\mathbf R\,\mathbf l\rvert^2=\lvert\mathbf S\mathbf R(\mathbf o-\mathbf c)\rvert^2\lvert\mathbf S\mathbf R\,\mathbf l\rvert^2-\left((\mathbf S\mathbf R(\mathbf o-\mathbf c))\cdot\mathbf S\mathbf R\,\mathbf l\right)^2$$ so $$\left(\lvert\mathbf S\mathbf R(\mathbf o-\mathbf c)\rvert^2-1\right) \lvert\mathbf S\mathbf R\,\mathbf l\rvert^2 \le \left((\mathbf S\mathbf R(\mathbf o-\mathbf c))\cdot\mathbf S\mathbf R\,\mathbf l\right)^2,$$ but $\mathbf Q=(\mathbf S\mathbf R)^T(\mathbf S\mathbf R)$, so this is just the original inequality derived from the quadratic equation, refactored.
Note that if you already have unit direction vectors for the principal axes, there’s no reason for the purposes of this calculation to decompose the rotation into three primitive rotations as in your question: $\mathbf R$ is simply the matrix with these vectors as its rows. Note also that since $\mathbf S$ is diagonal, you can implement multiplication of a vector by $\mathbf S$ with three multiplications.