As indicated in Erick's comment, your problem is not that the two calculations yield different results for a singular covariance matrix, but that $\mathbf S$ is singular (and hence not invertible) if some of the eigenvalues are zero, so that neither calculation is well-defined. This is a conceptual problem, not a computational one; the Mahalanobis distance is simply not well-defined in this case.
This paper suggests what it calls a regularized Mahalanobis distance to deal with this problem.
If anyone is still interested in this problem, in my spare time i think i found an answer. The answer is: the two methods only yield the same lengths for the semi-minor and semi-major-axis of the ellipse for the theoretical case that the covariance matrix is exactly equal to the expectation value of the covariance matrix for a bivariate gaussian distribution.
$$Cov = \left[{\begin{array}{cc}
\sigma_1^2 & \rho\sigma_1\sigma_2\\
\rho\sigma_1\sigma_2 & \sigma_2^2\\
\end{array} }\right]$$
If the covariance-matrix comes from a numerical realization of a bivariate gaussian distribution, then the covariance matrix will not exactly look like stated above and the results of the "eigenvalue"-method and the "Batschelet, 1981"-Method begin to deviate. The result can be seen in this python code below
import numpy
import scipy
# example 14.4.1 (p. 273), parameter for bivariate gaussian distribution
mu_1 = 5
mu_2 = 4
sigma_1 = 3.0
sigma_2 = 2.0
rho = 0.4 # correlation
# eq 14.4.3 (p. 269), expectation value for covariance-matrix for bivariate gaussian distribution!
cov_theo = numpy.array([[sigma_1**2, rho*sigma_1*sigma_2],
[rho*sigma_1*sigma_2, sigma_2**2]])
# make an example (with 10000 points) as sketched in Batschelet, 1981
points = numpy.random.multivariate_normal(mean=(mu_1,mu_2), cov=cov_theo, size=10000)
x, y = points.T
cov = numpy.cov(x, y) # covariance-matrix by the realization
# get estimations for mu, sigma, rho
mu_1 = numpy.mean(x)
mu_2 = numpy.mean(y)
sigma_1 = numpy.std(x)
sigma_2 = numpy.std(y)
rho = scipy.corrcoef(x,y)[0,1]
## -- BATSCHELET-METHOD (1981) -- ##
# eq 14.4. (ellipse equation): A(x-mu_1)^2 + 2B(x-mu_1)(y-mu_2) + C(y-mu_2)^2 = D
A = sigma_2**2
B = -rho*sigma_1*sigma_2
C = sigma_1**2
D = (1-rho**2)*sigma_1**2*sigma_2**2
# further calculations (p. 264)
R = numpy.sqrt((A-C)**2 + 4*B**2)
a = numpy.sqrt(2*D/(A+C-R))
b = numpy.sqrt(2*D/(A+C+R))
## -- EIGENVALUE-METHOD from SO -- ##
eigval, _ = numpy.linalg.eig(cov)
a_eig, b_eig = numpy.sqrt(eigval)
print 'realization of covariance_matrix\n', cov
print 'deviations from what it should look like:'
print 'sigma_1^2 - cov[0,0] = {d1}'.format(d1=sigma_1**2-cov[0,0])
print 'sigma_2^2 - cov[1,1] = {d2}'.format(d2=sigma_2**2-cov[1,1])
print 'rho*sigma_1*sigma_2 - cov[0,1] = {d}'.format(d=rho*sigma_1*sigma_2 - cov[0,1])
print 'a-a_eig = {diff_a}, b-b_eig = {diff_b}'.format(diff_a=a-a_eig, diff_b=b-b_eig)
The result looks like:
realization of covariance_matrix
[[ 9.00812576 2.4676689 ]
[ 2.4676689 4.14435948]]
deviations from what it should look like:
sigma_1^2 - cov[0,0] = -0.00090081257572
sigma_2^2 - cov[1,1] = -0.00041443594839
rho*sigma_1*sigma_2 - cov[0,1] = -0.000246766890107
a-a_eig = -0.000158440395525, b-b_eig = -8.82014644565e-05
and it gets worse for smaller sample sizes.
In Summary:
- The methods only give the same results, if the covariance matrix is exactly as assumed.
- The sample points must be drawn from a bivariate gaussian distribution. If the points are not taken from such a distribution, it gets even worse.
- I am still unable to show, that $a=\sqrt{\lambda_1}\quad,\quad b=\sqrt{\lambda_2}$ iff the covariance-matrix is exactly as assumed.
Best Answer
If I understand correctly your question, your formula seems to be wrong.
In three dimensions the distribution of mass could be represented by the inertia ellipsoid, whose equation, in a principal frame is $$ I_x x^2 + I_y y^2 + I_z z^2 = 1, $$ so that the length of the axes of the ellipsoid are $a_i = 1/\sqrt{I_i} = 1/\sqrt{\lambda_i}$ (the principal moments of inertia are the eigenvalues of the inertia matrix).
See section 4.5 of http://www.eng.auburn.edu/~marghitu/MECH2110/C_4.pdf .