I’m not an expert on MPS and MPO myself, but I am afraid this is not a trivial thing to do with computational efficiency.
However, for specific classes of mixed states that arise from partitions of a spin chain which is globally described by a matrix product state, this reference [1807.01640] Uhlmann fidelities from tensor networks presents some efficient methods.
I agree there probably isn’t an efficient way to compute this type of fidelity, given \rho and \sigma as MPOs. Is there a different type of fidelity that could work for the purpose you need?
I do know of one special set of cases where you can get the square root of an MPO, and that is where the MPO is the exponential of an operator H and you can make it by time evolving (real or imaginary) by H. Then you can just time evolve for half the time to get the square root.