Failing an analytic answer I would use Newton's method as a quick and dirty numerical iterator:
You can rearrange to:
$$ (1-\rho) \left (\Phi^{-1}(LGD.DR) - \Phi^{-1}(DR) \right ) + \Phi^{-1}(PD) - \Phi^{-1}(PD.LGD) = 0$$
which given you have fixed DR, LGD and $\rho$ is essentially:
$$K + \Phi^{-1}(PD) - \Phi^{-1}(PD.LGD) = f(PD) = 0$$
Newton's iterative formula yeilds the iterative scheme:
$$ PD_{n+1} = PD_n + \frac{f(PD_n)}{f'(PD_n)} $$
or $$PD_{n+1} = PD_n + \frac{K + \Phi^{-1}(PD_n) - \Phi^{-1}(PD_n.LGD)}{\frac{\partial}{\partial PD} \left( \Phi^{-1}(PD_n) - \Phi^{-1}(PD_n.LGD) \right )} $$
And actually there are numerous algorithms you might use - its not a computationally intensive task: in python you can try:
https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.root_scalar.html#scipy.optimize.root_scalar
from scipy import optimize
def f(x, *args):
val = # define your function in terms of x = PD, and other static args
return val
sol = optimize.root_scalar(f, args={}, bracket=[0, 3], method='brentq')
sol = optimize.root_scalar(f, args={}, x0=0.2, fprime=fprime, method='newton')