Let's say $X \sim N(0, 1)$ is a standard Gaussian input signal and $Y \in \{y_1, y_2, \cdots, y_n\}$ is a discrete reconstruction alphabet. Under the MSE distortion measure $d(x, y) = (x - y)^2$, for a given rate $R^*$, I want to find

$$ \begin{aligned} D(R^*) =& \inf_{p(y|x)} \int p(x) p(y|x) d(x, y) dxdy \\ & \text{subject to}\; I(X; Y) \le R^* \end{aligned}. $$

This is a convex optimization problem which is generally easy, but here $p(y|x)$ can be an arbitrary continuous distribution so I don't know how to parameterize it. I am aware of the Blahut–Arimoto algorithm, but it only works when both $X$ and $Y$ are discrete and finite (in which case $p(y|x)$ can be expressed as a table).

Anyway, it sounds like a pretty fundamental question. Is there a known result about it?