The problem of reducing a fully anisotropic (triclinic) stiffness tensor comprised of 21 distinct components to a higher symmetry form having fewer distinct elements is of interest in many geophysical applications, where assuming a particular type of symmetry is often necessary to sufficiently reduce computational complexity to allow practical solutions. In addition, recent advances in the upscaling of realistically large field-scale discrete fracture networks have led to the need for an efficient way to derive a very large number of nearest medium approximations to the triclinic stiffness tensors obtained. Owing to rotational symmetries and nonlinearity, the problem of efficiently finding such approximations is generally nontrivial because optimal orientations are intrinsically nonunique. An algorithm proposed by Dellinger computes nearest orthotropic and transverse isotropic approximations for a given stiffness tensor, using the Federov norm as an objective function to iteratively minimize the fit error. Although this method is appropriate for computing solutions to single problem instances, the implementation is too inefficient for production situations, where a very large number of invocations of the algorithm is required. The enhanced algorithm proposed here is accurate, efficient, and general, allowing nearest medium approximations to be determined for arbitrary symmetry types, including isotropic, cubic, transverse isotropic, orthotropic, and monoclinic.