I think a more fundamental way to approach the problem is by discussing geodesic curves on the surface you call home. Remember that the geodesic equation, while equivalent to the Euler-Lagrange equation, can be derived simply by considering differentials, not extremes of integrals. The geodesic equation emerges exactly by finding the acceleration, and hence force by Newton's laws, in generalized coordinates.
See the Schaum's guide Lagrangian Dynamics by Dare A. Wells Ch. 3, or Vector and Tensor Analysis by Borisenko and Tarapov problem 10 on P. 181
So, by setting the force equal to zero, one finds that the path is the solution to the geodesic equation. So, if we define a straight line to be the one that a particle takes when no forces are on it, or better yet that an object with no forces on it takes the quickest, and hence shortest route between two points, then walla, the shortest distance between two points is the geodesic; in Euclidean space, a straight line as we know it.
In fact, on P. 51 Borisenko and Tarapov show that if the force is everywhere tangent to the curve of travel, then the particle will travel in a straight line as well. Again, even if there is a force on it, as long as the force does not have a component perpendicular to the path, a particle will travel in a straight line between two points.
Also, as far as intuition goes, this is also the path of least work.
So, if you agree with the definition of a derivative in a given metric, then you can find the geodesic curves between points. If you define derivatives differently, and hence coordinate transformations differently, then it's a whole other story.
Best Answer
The surface is $$z(x, y) = a x^2 + b y^2 + c x y + d x + e y + f \tag{1}\label{AC1}$$ The axis of symmetry is parallel to the $z$ axis, and $z(x, y)$ has a stationary point where $(x, y)$ are at the axis of symmetry: $$\left\lbrace \begin{aligned} \displaystyle \frac{\partial z(x, y)}{\partial x} &= 0 \\ \displaystyle \frac{\partial z(x, y)}{\partial y} &= 0 \\ \end{aligned} \right. $$ This has exactly one solution, $$\left\lbrace \begin{aligned} \displaystyle x_0 &= \frac{2 b d - c e}{c^2 - 4 a b} \\ \displaystyle y_0 &= \frac{2 a e - c d}{c^2 - 4 a b} \\ \end{aligned}\right. \tag{2a}\label{AC2a}$$ Note that $$z_0 = z(x_0 , y_0) = f + \frac{a e^2 - c d e + b d^2}{c^2 - 4 a b} \tag{2b}\label{AC2b}$$ and that if $c^2 = 4 a b$, the surface is not a paraboloid, as it does not have an axis of symmetry.
To make the next steps clearer to follow, let's reparametrize the surface as $(u, v)$: $$\left\lbrace\begin{aligned} x(u, v) &= u \\ y(u, v) &= v \\ z(u, v) &= a u^2 + b v^2 + c u v + d u + e v + f \\ \end{aligned}\right. \tag{3}\label{AC3}$$ The distance squared between point $(x_i , y_i , z_i)$ and $(u, v) = (x(u,v), y(u,v), z(u,v))$ is $s(u, v)$, $$s(u, v) = (u - x_i)^2 + (v - y_i)^2 + (a u^2 + b v^2 + c u v + d u + e v + f - z_i)^2 \tag{4a}\label{AC4a}$$ i.e. $$\begin{array}{rl} s(x, y) & = (a^2) x^4 + (b^2) y^4 \\ ~ & + (2 a c) x^3 y + (2 b c) x y^3 \\ ~ & + (2 a d) x^3 + (2 b e) y^3 \\ ~ & + (c^2 + 2 a b) x^2 y^2 \\ ~ & + (2 a e + 2 c d) x^2 y \\ ~ & + (2 c e + 2 b d) x y^2 \\ ~ & + (1 + d^2 + 2 a (f - z_i)) x^2 \\ ~ & + (1 + e^2 + 2 b (f - z_i)) y^2 \\ ~ & + (2 d e + 2 c (f - z_i)) x y \\ ~ & + (2 d (f - z_i) - 2 x_i) x \\ ~ & + (2 e (f - z_i) - 2 y_i) y \\ ~ & + (x_i^2 + y_i^2 + (f - z_i)^2) \\ \end{array} \tag{4b}\label{AC4b}$$ Its partial derivatives define the gradient, $$\left\lbrace ~ \begin{array}{rl} \displaystyle \frac{\partial s(u,v)}{\partial u} = s_u(u,v) & = (4 a^2) u^3 ~ + (2 b c) v^3 \\ ~ & + (6 a c) u^2 v ~ + (2 c^2 + 4 a c) u v^2 \\ ~ & + (6 a d) u^2 ~ + (2 c e + 2 b d) v^2 \\ ~ & + (4 a e + 4 c d) u v \\ ~ & + (2 + 2 d^2 + 4 a f - 4 a z_i) u \\ ~ & + (2 d e + 2 c f - 2 c z_i) v \\ ~ & + (2 d f - 2 d z_i - 2 x_i) \\ \displaystyle \frac{\partial s(u,v)}{\partial v} = s_v(u,v) & = (2 a c) u^3 ~ + (4 b^2) v^3 \\ ~ & + (2 c^2 + 4 a b) u^2 v ~ + (6 b c) u v^2 \\ ~ & + (2 a e + 2 c d) u^2 ~ + (6 b e) v^2 \\ ~ & + (4 c e + 4 b d) u v \\ ~ & + (2 c f + 2 d e - 2 c z_i) u \\ ~ & + (2 + 2 e^2 + 4 b f - 4 b z_i) v \\ ~ & + (2 e f - 2 e z_i - 2 y_i) \\ \end{array}\right. \tag{5}\label{5}$$ Although the distance squared function has 15 coefficients (up to for $u^4$ and for $v^4$), it is a surprisingly smooth function. Starting at $u = x_i$ and $v = y_i$, gradient descent should find the true minimum.
If we call the point on the surface $(u_i, v_i) = (u + u_{\Delta}, v + v_{\Delta})$ closest to the point $(x_i , y_i, z_i)$, we have immediately limited the region of interest to $u_{\Delta}^2 + v_{\Delta}^2 \le s(u, v)$, because all other points on the surface are more distant from surface point $(u, v)$ than the point of interest $(x_i , y_i , z_i)$.
Finding the minimum distances for a million points $i$ should not take that long on current computers; a few seconds, I'd guess.
One way to speed up the calculation is to choose better starting points $(u, v)$ for the iterations (instead of simply $u = x_i$, $v = y_i$). A reasonable method would be to set up a regular rectangular grid, and find the $u, v$ coordinates for the closest point on the surface to the center of that cell. (You would still calculate the squared distance to both points on the surface, and start the gradient descent or other iterative method at the one with the smaller squared distance.) Note that you can easily set up a sentinel coordinate pair to indicate "unknown", and only find the coordinates when needed; this is important, if the points are clumped within the grid, with many grid cells without any points at all. Still, this method is only useful if each grid cell has many points inside it, on average.
When doing the gradient descent, remember that all but $x$ and $y$ are constants, and you can precalculate the coefficients to speed up the evaluation at each candidate point. (To evaluate the squared distance, and its partial derivatives, you'll need around 44 multiplications and around 32 additions or subtractions.)
It is possible to translate and rotate the coordinate system (around the $z$ axis, in the $x y$ plane) so that the surface is just $$z(x, y) = C_{20} x^2 + C_{02} y^2$$ and the squared distance between point $z(x, y)$ on the surface and point $(x_i , y_i , z_i)$ is $$s(x, y) = (x - x_i)^2 + (y - y_i)^2 + (C_{20} x^2 + C_{02} y^2 - z_i)^2$$ and its partial derivatives are $$\begin{aligned} s_x (x,y) &= (4 C_{20}^2) x^3 + (4 C_{20} C_{02}) x y^2 + (2 - 4 C_{20} z_i) x - 2 x_i \\ s_y (x,y) &= (4 C_{02}^2) y^3 + (4 C_{20} C_{02}) x^2 y + (2 - 4 C_{02} z_i) y - 2 y_i \\ \end{aligned}$$ For each candidate point $(x, y)$, you'll need to do seven multiplications and six additions or subtractions to get the squared distance; or a total of seventeen multiplications and fourteen additions or subtractions for the squared distance and the partial derivatives. This can be a significant speedup, if the gradient descent requires many steps to find the minimum.
If we use a translation by $(x_0, y_0, z_0)$ and a rotation of $\theta$ counterclockwise on the $x y$ plane, the original point $(x_i , y_i , z_i)$ is transformed to $$\left\lbrace ~ \begin{aligned} x_k &= x_0 + x_i \cos\theta - y_i \sin\theta \\ y_k &= y_0 + y_i \cos\theta + y_i \sin\theta \\ z_k &= z_0 + z_i \\ \end{aligned} \right. \iff \left\lbrace ~ \begin{aligned} x_i &= x_k \cos\theta + y_k \sin\theta - x_0 \\ y_i &= y_k \cos\theta - x_k \sin\theta - y_0 \\ z_i &= z_k - z_0 \\ \end{aligned} \right. $$
The surface in the transformed coordinates (subscript $k$) is $$z = C_{20} x^2 + C_{02} y_2$$ which in original coordinates corresponds to $$\begin{aligned} z + z_0 & = C_{20} \left( x_0 + x_i \cos\theta - y_i \sin\theta \right)^2 \\ & + C_{02} \left( y_0 + y_i \cos\theta + x_i \sin\theta \right)^2 \\ \end{aligned}$$ ie. $$\begin{aligned} z(x, y) & = ( C_{02} (\sin\theta)^2 + C_{20} (\cos\theta)^2 ) x^2 \\ ~ & + ( C_{20} (\sin\theta)^2 + C_{02} (\cos\theta)^2 ) y^2 \\ ~ & + ( 2 (C_{02} - C_{20}) \cos\theta\sin\theta ) x y \\ ~ & + ( 2 C_{20} x_0 \cos\theta - 2 C_{02} y_0 \sin\theta ) x \\ ~ & + ( 2 C_{02} y_0 \cos\theta - 2 C_{20} x_0 \sin\theta ) y \\ ~ & + ( C_{20} x_0^2 + C_{02} y_0^2 ) \\ \end{aligned}$$ Comparing the coefficients to the original form, we see that $$\begin{aligned} a &= C_{02} (\sin\theta)^2 + C_{20} (\cos\theta)^2 \\ b &= C_{20} (\sin\theta)^2 + C_{02} (\cos\theta)^2 \\ c &= ( 2 C_{02} - C_{20}) \cos\theta\sin\theta \\ d &= 2 C_{20} x_0 \cos\theta - 2 C_{02} y_0 \sin\theta \\ e &= 2 C_{02} y_0 \cos\theta - 2 C_{20} x_0 \sin\theta \\ f &= C_{20} x_0^2 + C_{02} y_0^2 \\ \end{aligned}$$