Sanity check: the statement is indeed true when $\rho$ is a pure state.
We can start by finding the singular values of the combination of purified systems, which I will write as $|\phi\rangle$ and $|\Psi\rangle$. Given the Hermitian matrix $M=|\phi\rangle\langle\phi|-|\Psi\rangle\langle\Psi|$, we can look for eigenvectors of the form $\alpha|\phi\rangle+\beta|\Psi\rangle$ (this is much easier when we have a sum of projectors, not a difference like we have here). Defining $z=\langle\phi|\Psi\rangle$, the eigenvalue equation reads
\begin{align}
M\left(\alpha|\phi\rangle+\beta|\Psi\rangle\right)&=|\phi\rangle\left(\alpha+\beta z\right)+|\Psi\rangle\left(-\beta-\alpha z^*\right)\\
&=\lambda\left(\alpha|\phi\rangle+\beta|\Psi\rangle\right).
\end{align} The eigenvalues are given by
$$\lambda=\frac{\alpha+\beta z}{\alpha}=-\frac{\beta+\alpha z^*}{\beta},$$ which means that the eigenvectors must obey
$$\alpha=\beta\frac{-1\pm\sqrt{1-|z|^2}}{z^*};$$ we thus have the unnormalized eigenvectors and eigenvalues
$$|v_\pm\rangle=\left(-1\pm\sqrt{1-|z|^2}\right)|\phi\rangle+z^*|\Psi\rangle,\qquad \lambda_\pm=\mp\sqrt{1-|z|^2}.$$ The entire matrix $M$ thus has two degenerate singular values $$\sigma_1=\sigma_2=\sqrt{1-\left|\langle\phi|\Psi\rangle\right|^2}.$$ This means that we can write your purification's trace distance as $$\big|\big||\phi\rangle\langle\phi|-|\Psi\rangle\langle\Psi|\big|\big|_1=2\sqrt{1-\left|\langle\phi|\Psi\rangle\right|^2}.$$ To achieve the smallest value of $\big|\big||\phi\rangle\langle\phi|-|\Psi\rangle\langle\Psi|\big|\big|_1$, we require the largest value of $\left|\langle\phi|\Psi\rangle\right|^2$.
Now, this expression for the singular values uses the overlap between two purified states, so we can connect it to the fidelity! Labelling all viable purifications of $\rho$ by $|\phi_\rho\rangle$, we know that $$F(\rho,\psi)=\max_{|\phi_\rho\rangle}\left|\langle\phi_\rho|\Psi\rangle\right|^2,$$ where we have fixed our purification $|\Psi\rangle$ of $|\psi\rangle$ without loss of generality. So we can collect our results and write $$\min_{|\phi_\rho\rangle}\big|\big||\phi_\rho\rangle\langle\phi_\rho|-|\Psi\rangle\langle\Psi|\big|\big|_1=2\sqrt{1-F(\rho,\psi)}\stackrel{?}{\leq}\epsilon.$$
Our question is now how to relate $F(\rho,\psi)$ and $\epsilon$. We know that $1-F(\rho,\psi)\leq {\epsilon}/2$, so the question is by how much can it be smaller. We can always write this fidelity as $F(\rho,\psi)=\langle \psi|\rho|\psi\rangle.$ Choosing some fiducial state $\rho=p|\psi\rangle\langle\psi|+(1-p)|\psi_\perp\rangle\langle\psi_\perp|$ with orthonormal $|\psi\rangle$ and $|\psi_\perp\rangle$, we can calculate (2 times) the trace distance $$\big|\big|\rho-|\psi\rangle\langle\psi|\big|\big|_1=2\left(1-p\right)$$ and the fidelity $$F(\rho,\psi)=p.$$ Since $$2\sqrt{1-p}\geq 2\left(1-p\right),$$ we have a counterexample where $$\min_{|\phi_\rho\rangle}\big|\big||\phi_\rho\rangle\langle\phi_\rho|-|\Psi\rangle\langle\Psi|\big|\big|_1 > \big|\big|\rho-|\psi\rangle\langle\psi|\big|\big|_1,$$ so it does not seem like the desired property always holds for all $\epsilon$. In fact, it seems as though this counterexample saturates the $\sqrt{\epsilon}$ condition that we have been trying to beat.
