Mathematically, the quantity you're asking about is called the operator norm. Unfortunately, there's no simple formula for it. If it's a fully general affine transformation - for instance, if it could have an arbitrary combination of rotations and nonuniform scales, in any order - then I'm afraid there's nothing for it but to use singular value decomposition. If you apply SVD to your matrix then the largest singular value will be the maximum radius of the resulting ellipsoid. The other singular values will also be its other two radii, and the SVD procedure can also extract the orientation of the axes for you.
If the form of your matrix is restricted in such a way that nonuniform scale happens at most once and is the first transform applied, i.e. is rightmost when you're using column vectors, then you can simplify this to just look at the lengths of the transformed axis vectors. In that case alone - i.e. a single nonuniform scale followed by any sequence of rotations, reflections, and uniform scales - looking at just the axis vectors will give you the right answer.