inverse bilinear map is a good example where, even if there is an analytical solution, it's often more robust to just do a couple newton iterations until you get down to machine precision.
added benefit is that it carries on in 3D (although with a larger jacobian of course).
You are absolutely right re: robustness (and that code falls back to the iterative method), but IIRC the Newton scheme was at least one order of magnitude slower than the formula. It was also slightly less accurate simply because it did more floating point operations.
I have not yet figured out how to make the 3D formula work so we always do the iterative method in 3D.
A nice trick: we approximate the (bi/tri)linear map as affine (usually quite accurate) to get a good initial guess.
added benefit is that it carries on in 3D (although with a larger jacobian of course).