Edit: my original solution didn't work. The eigenvalues of the mappping matrix only correspond to the scale factors for special orientation cases.
Fortunately, I've found a solution that seems to hold in general.
If we take a collection of unit vectors in UV space and map them into 3D, they will form an ellipse, whose major and minor axes correspond to the stretch factors we want.
We don't care about the absolute orientation of these vectors after transformation, just how their length changes, according to the equation of an ellipse centered at the origin:
A * u*u + B * u*v + C * v*v = l*l
So first we need to find the parameters A, B, and C. (I'll rename the vertices to V1, V2, V3 so they don't conflict with the conventional labeling of these coefficients)
Consider the set of edges of the triangle (in both 3D and uv space):
E1 = V1 - V3
E2 = V2 - V1
E3 = V3 - V2
We can construct a system of equations to find the parameters as follows:
E1.u * E1.u & E1.u * E1.v & E1.v * E1.v \\
E2.u * E2.u & E2.u * E2.v & E2.v * E2.v \\
E3.u * E3.u & E3.u * E3.v & E3.v * E3.v \\
E1.xyz\cdot E1.xyz \\
E2.xyz\cdot E2.xyz \\
E3.xyz\cdot E3.xyz \\
Note that by construction this is insensitive to the rotation of the triangle in 3D, which is a good start, and doesn't need any normalizations to get there unlike my previous solution.
If the determinant of the coefficient matrix is zero, then the UV triangle is degenerate (at least one of the scale factors is infinite), and needs special-case handling which I'll elide for now.
After solving this by your favourite method (I used Cramer's Rule for a crude test), we have two cases for finding the axis scale factors:
// Ellipse is symmetrical across the line u=v.
// (Either theta = PI/4, or B = 0 and scaling is uniform).
p = A - 0.5f * B;
q = A + 0.5f * B;
// Ellipse is non-uniformly scaled and arbitrarily rotated.
// (Or a rounding error makes it appear slightly so)
// Find its angle:
float theta = atan(B/(C - A))/2f;
// cos squared
float c = cos(theta);
c *= c;
// sin squared
float s = 1f - c;
float divisor = 1f/(c - s);
p = (A * c - C * s) * divisor;
q = (C * c - A * s) * divisor;
// Note that numerical errors can sometimes cause p & q to be slightly negative.
// This isn't meaningful, so clamp to a non-negative value before taking the root.
float scale1 = sqrt(pos(p));
float scale2 = sqrt(pos(q));
Although I didn't need it, this calculates the direction of non-uniform stretch too (theta). It does not directly detect reflection, but this can be added by checking the winding of the triangle in uv space.