Fix $(\varphi, \psi) \in \Phi_c$.
- There are a $\mu$-null set $N_x$ and $\nu$-null net $N_y$ such that $\varphi(x)+\psi(y) \le c(x,y)$ for all $(x,y) \in N_x^c \times N_y^c$.
- We re-define $\varphi, \psi$ by letting them take value $-\infty$ on $N_x, N_y$ respectively. In this way, $\varphi(x)+\psi(y) \le c(x,y)$ for all $(x,y) \in X \times Y$.
We define $\varphi^c$ by
$$
\varphi^c (y) := \inf_{x\in X} [c(x,y) - \varphi(x)].
$$
- The infimum of a collection of extended real-valued measurable functions is again a measurable function, so $\varphi^c$ is measurable.
- There is $x_0 \in N^c_x$, so $\varphi^c (y) \le c(x_0,y)-\varphi(x_0) <+\infty$ for all $y \in Y$.
- Given $y\in Y$, $\varphi(x)+\psi(y) \le c(x,y)$ for all $x \in X$, so $\varphi^c (y) \ge \psi(y)$ for all $y\in Y$.
We define $\varphi^{cc}$ by
$$
\varphi^{cc} (x) := \inf_{y\in Y} [c(x,y) - \varphi^c (y)].
$$
For all $x\in X$,
$$
\begin{align}
\varphi^{cc} (x) &= \inf_{y\in Y} \left [c(x,y) - \inf_{z\in X} [c(z,y) - \varphi(z)] \right] \\
&= \inf_{y\in Y} \left [c(x,y) + \sup_{z\in X} [-c(z,y) + \varphi(z)] \right] \\
&= \inf_{y\in Y} \sup_{z\in X} [c(x,y) -c(z,y) + \varphi(z)] \\
&\ge \inf_{y\in Y} \varphi(x) =\varphi(x) \quad \text{by picking} \quad z=x.
\end{align}
$$
There is $y_0 \in N_y^c$. Then
$$
\varphi^{cc} (x) \le c(x, y_0) - \varphi^c (y_0) \le c(x, y_0) - \psi(y_0) <+\infty \quad \forall x\in X.
$$
For all $(x,y) \in X \times Y$,
$$
\varphi^{cc} (x) + \varphi^c (y) = \inf_{z\in Y} [c(x,z) - \varphi^c(z)] + \varphi^c (y) \le [c(x,y) - \varphi^c (y)] + \varphi^c (y)= c(x,y).
$$
- There are a $\mu$-null set $M_x$ and $\nu$-null set $M_y$ such that $c (x, y) \le c_X(x)+c_Y(y)$ for all $(x, y) \in M_x^c \times M_y^c$.
- We re-define $c_X, c_Y$ such that $c_X (x) = c_Y(y) := +\infty$ for all $(x, y) \in M_x \times M_y$. In this way, $c(x,y) \le c_X(x)+c_Y(y)$ for all $(x,y) \in X \times Y$.
Let
$$
a := \inf_{y\in Y} [c_Y(y) - \varphi^c (y)].
$$
First,
$$
a \le c_Y(y_0) - \varphi^c (y_0) \le c_Y(y_0) - \psi (y_0) < +\infty.
$$
There is $x_1 \in (N_x \cup M_x)^c$, so
$$
\begin{align}
c_Y(y) - \varphi^c (y) &= c_Y(y) - \inf_{x\in X} [c(x,y) - \varphi(x)] \\
&= \sup_{x\in X} [c_Y(y)-c(x,y) + \varphi(x)] \\
&\ge \sup_{x\in X} [\varphi(x)-c_X(x)] \\
&\ge \varphi(x_1)-c_X(x_1) \\
&> -\infty.
\end{align}
$$
It follows that $a \in \mathbb R$. Clearly, $a \le c_Y(y) - \varphi^c (y)$ for all $y\in Y$, so $\varphi^c+a \le c_Y$. We have
$$
\begin{align}
(\varphi^{cc} (x)-a) - c_X(x) &= \inf_{y\in Y} [c(x,y) - \varphi^c (y)]-a - c_X(x) \\
&= \inf_{y\in Y} [(c(x,y)-c_X(x)) - \varphi^c (y)] -a \\
&\le \inf_{y\in Y} [c_Y(y) - \varphi^c (y)] -a \\&=0.
\end{align}
$$
Clearly, $(\varphi', \psi') := (\varphi^{cc}-a, \varphi^c+a)$ satisfies the requirement.
Let $c_X = c_Y := \frac{1}{2} |\cdot|^2$. We need the following Lemma
Lemma: If $X,Y$ are Polish spaces, $c$ l.s.c., then
$$
\min_{\pi \in \Pi(\mu, \nu)} \mathbb K (\pi) = \sup _{(\varphi, \psi) \in\Phi_c} \mathbb J (\varphi, \psi) .
$$
If, moreover, there are $(c_X, c_Y) \in L_1(\mu) \times \in L_1 (\nu)$ such that $c (x, y) \le c_X(x)+c_Y(y)$ for $\mu$-a.e. $x\in X$ and $\nu$-a.e. $y\in Y$, then the maximization
$$
\sup_{(\varphi, \psi) \in \Phi_c} J (\varphi, \psi)
$$
has a solution of the form $(\psi, \psi^c)$ for some $\psi \in L_1(\mu)$ with $\psi^{cc} = \psi$. Here $\psi^c$ is the $c$-transform of $\psi$.
Let $\pi^\dagger$ minimize $\mathbb K$ over $\Pi(\mu, \nu)$. By Lemma, there is $(\varphi, \varphi^c)$ with $\varphi^{cc} = \varphi$ maximizing $\mathbb J$ over $\Phi_c$. Because $c(x, y) = c_X(x) + c_Y(y) - \langle x, y\rangle$, we have $(\tilde \varphi, \tilde \psi) := (c_X - \varphi, c_Y - \varphi^c)$ minimizes $\mathbb J$ over $\tilde \Phi$. It's can be shown that $\tilde \psi$ is the convex conjugate of $\tilde \varphi$, i.e., $\tilde \psi = \tilde \varphi^*$, and that $\tilde \varphi$ is the convex conjugate of $\tilde \varphi^*$, i.e., $\tilde \varphi = \tilde \varphi^{**}$. By this result, $\tilde \varphi$ is convex l.s.c. By Lemma,
$$
\int_X \tilde \varphi (x) d \mu (x) + \int_Y \tilde \psi (y) d \nu (y)= \int_{X \times Y} \langle x, y\rangle d\pi^\dagger(x, y),
$$
or equivalently
$$
\int_{X \times Y} [\tilde \varphi (x) + \tilde \varphi^* (y) - \langle x, y\rangle] d\pi^\dagger(x, y) = 0.
$$
We have $\tilde \varphi (x) + \tilde \varphi^* (y) - \langle x, y\rangle \ge 0$ for all $(x,y) \in X \times Y$, so $\tilde \varphi (x) + \tilde \varphi^* (y) - \langle x, y\rangle = 0$ for $\pi^\dagger$-a.e. $(x, y) \in X \times Y$. By this result, $y \in \partial \tilde \varphi (x)$ for $\pi^\dagger$-a.e. $(x, y) \in X \times Y$. It follows from $\tilde \varphi = \tilde \varphi^{**}$ that $\tilde \varphi$ is convex l.s.c.
Let $\pi^\dagger \in \Pi(\mu, \nu)$ and $\tilde \varphi \in L_1 (\mu)$ be such that $y \in \partial \tilde \varphi (x)$ for $\pi^\dagger$-a.e. $(x, y) \in X \times Y$. Then $\tilde \varphi (x) + \tilde \varphi^* (y) - \langle x, y\rangle = 0$ for $\pi^\dagger$-a.e. $(x, y) \in X \times Y$, so
$$
\int_{X \times Y} [\tilde \varphi (x) + \tilde \varphi^* (y) - \langle x, y\rangle] d\pi^\dagger(x, y) = 0.
$$
By this result, $\langle \cdot, \cdot \rangle$ is $\pi^\dagger$-integrable, so the map $(x, y) \mapsto \tilde \varphi (x) + \tilde \varphi^* (y)$ is also $\pi^\dagger$-integrable. Because $\tilde \varphi$ is $\mu$-integrable, so $\tilde \varphi^* \in L_1 (\nu)$. It follows that
$$
\int_X\tilde \varphi d \mu + \int_Y \tilde \varphi^* d \nu = \int_{X \times Y}\langle x, y\rangle d\pi^\dagger(x, y).
$$
On the other hand,
$$
\int_X \varphi' d \mu + \int_Y \psi' d \nu \ge \int_{X \times Y} \langle x, y\rangle d\pi^\dagger(x, y) \quad \forall (\varphi', \psi') \in \tilde \Phi.
$$
It follows that $(\tilde \varphi, \tilde \varphi^*)$ minimizes $\mathbb J$ over $\tilde \Phi$. Because $c(x, y) = c_X(x) + c_Y(y) - \langle x, y\rangle$, we have $(\varphi, \psi) := (c_X - \tilde\varphi, c_Y - \tilde \varphi^*)$ minimizes $\mathbb J$ over $\Phi_c$. This completes the proof.
Best Answer
It's clear that $$ \sup _{\Phi_{c} \cap C_{b}} J(\varphi, \psi) \leq \sup _{\Phi_{c} \cap L^{1}} J(\varphi, \psi) \leq \inf _{\Pi(\mu, \nu)} I[\pi]. $$
So it remains to prove $$ \inf _{\Pi(\mu, \nu)} I[\pi] \le \sup _{\Phi_{c} \cap C_{b}} J(\varphi, \psi) . $$
To simplify notations, let $\varphi \oplus \psi: (x,y) \mapsto \varphi(x)+\psi(y)$. We have $$ \inf _{\pi \in \Pi(\mu, \nu)} I[\pi]=\inf _{\pi \in M_{+}(X \times Y)}\left(I[\pi]+\left\{\begin{array}{l} 0 \text { if } \pi \in \Pi(\mu, \nu) \\ +\infty \text { else } \end{array}\right)\right. $$ with $M_+(X \times Y)$ the space of non-negative Borel measures on $X\times Y$. Also, $$ \left\{\begin{array}{l} 0 \text { if } \pi \in \Pi(\mu, \nu) \\ +\infty \text { else } \end{array}\right\}=\sup \left \{ \int \varphi d \mu+\int \psi d \nu-\int \varphi \oplus \psi d \pi \,\middle\vert\, (\varphi, \psi) \in C_b(X) \times C_b(Y) \right \}, $$
It follows that $$ \begin{aligned} &\inf _{\pi \in \Pi(\mu, \nu)} I[\pi]\\ =&\inf _{\pi \in M_{+}(X \times Y)} \sup \left \{ \int_{X \times Y} c d \pi +\int_{X} \varphi d \mu+\int_{Y} \psi d \nu - \int_{X \times Y} \varphi \oplus \psi d \pi \,\middle\vert\, (\varphi, \psi) \in C_b(X) \times C_b(Y) \right\}. \end{aligned} $$
Let $E:=C_{b}(X \times Y)$ be the set of all bounded continuous functions on $X \times Y$, equipped with its usual supremum norm $\|\cdot\|_{\infty}$.
Because $X,Y$ are compact, so is $X\times Y$. So $E$ coincides with the space $C_{0}(X \times Y)$ of all continuous functions vanishing at infinity on $X \times Y$.
By Riesz' theorem, the topological dual $E^*$ of $E$ can be identified with the space of regular Radon measures, $M(X \times Y)$, normed by total variation.
Moreover, a nonnegative linear form on $E$ corresponds with a regular nonnegative Borel measure.
Then we introduce $$ \begin{gathered} \Theta: u \in E \longmapsto\left\{\begin{array}{l} 0 &\text {if } u \geq-c, \\ +\infty &\text {else}. \end{array}\right.\\ \text{}\\ \Xi: u \in E \longmapsto\left\{\begin{array}{l} \int_{X} \varphi d \mu+\int_{Y} \psi d \nu &\begin{align*} &\text {if } u = \varphi \oplus \psi \text{ for}\\ & \text{some } (\varphi, \psi) \in C_b(X) \times C_b(Y) \end{align*},\\ +\infty &\text {else. } \end{array}\right. \end{gathered} $$
It's easy to verify that $\Theta, \Xi$ satisfy the condition of Fenchel-Rockafellar duality, so $$ \inf _{u\in E}[\Theta(u)+\Xi(u)] = \max _{\pi \in E^{*}}\left[-\Theta^{*}\left(-\pi\right)-\Xi^{*}\left(\pi\right)\right]. $$
It's clear that $$ \begin{align*} \inf _{u\in E}[\Theta(u)+\Xi(u)] &= \inf \left\{\int_{X} \varphi d \mu+\int_{Y} \psi d \nu \,\middle\vert\, (\varphi, \psi) \in C_b(X) \times C_b(Y) \text{ s.t. } \varphi \oplus \psi \geq -c \right\} \\ &=-\sup \left\{J(\varphi, \psi) \mid (\varphi, \psi) \in \Phi_{c}\cap C_{b}\right\}. \end{align*} $$
Next, we compute the Legendre-Fenchel transforms of $\Theta, \Xi$. First, for any $\pi \in E^*$, $$ \begin{aligned} \Theta^{*}(-\pi) =\sup _{u \in E}\left\{-\int u d \pi \,\middle\vert\, u \geq-c\right\} = \sup _{u \in E}\left\{\int u d \pi \,\middle\vert\, u \leq c\right\} . \end{aligned} $$
Thus $$ \Theta^{*}(-\pi) = \begin{cases} \int c d \pi &\text {if } \pi \in M_{+}(X \times Y) \\ +\infty &\text {else}. \end{cases} $$
We also have $$ \begin{align*} \Xi^{*}(\pi) &= \sup_{u\in E} \left \{ \int ud\pi - \int \varphi d \mu- \int \psi d\nu \,\middle\vert\, u = \varphi \oplus \psi \text{ for some } (\varphi, \psi) \in C_b(X) \times C_b(Y) \right \} \\ &= \begin{cases} 0 &\text {if } \quad \forall(\varphi, \psi) \in C_b(X) \times C_b(Y) : \int \varphi \oplus \psi d \pi= \int \varphi d \mu+\int \psi d \nu \\ +\infty & \text {else} \end{cases} \\ &= \begin{cases} 0 &\text {if } \quad \pi \in \Pi(\mu, \nu) \\ +\infty & \text {else}. \end{cases} \end{align*} $$
It follows that $$ \max _{\pi \in E^{*}}\left[-\Theta^{*}\left(-\pi\right)-\Xi^{*}\left(\pi\right)\right] = \max _{\pi \in \Pi(\mu, \nu) \cap M_{+}(X \times Y)} -\int c d \pi = - \min _{\pi \in \Pi(\mu, \nu)} \int c d \pi. $$ Hence $$ -\sup \left\{J(\varphi, \psi) \mid (\varphi, \psi) \in \Phi_{c} \cap C_{b}\right\} = - \min _{\pi \in \Pi(\mu, \nu)} \int c d \pi. $$
Fix $$ \pi_* \in \operatorname{argmin}_{\pi \in \Pi(\mu, \nu)}I [\pi]. $$
Then $\pi_*$ is tight, so there are compact sets $X_n \subset X$ and $Y_n \subset Y$ such that $\mu(X_n^c), \nu(Y_n^c) \le \delta$. Let $Z_n := X_n \times Y_n$. Then $\pi (Z_n^c) \le 2\delta$. We define $\pi_n \in \mathcal P(Z_n)$ by $$ \pi_n := \frac{1_{Z_n}\pi_*}{\pi_* (Z_n)}. $$
Let $\mu_n := P_\sharp^{X_n} \pi_{n} \in \mathcal P(X_n)$ and $\nu_n := P_\sharp^{Y_n} \pi_{n} \in \mathcal P(Y_n)$ be marginals of $\pi_n$. Let $$ I_n[\pi] := \int_{Z_n} c d \pi \quad \forall \pi \in \Pi(\mu_n, \nu_n). $$
Fix $$ \tilde \pi_{n} \in \operatorname{argmin}_{\pi \in \Pi(\mu_n, \nu_n)} I_n [\pi]. $$
Define $\pi_{*n}$ by $$ \pi_{*n} := \pi_* (Z_n) \tilde \pi_{n} + 1_{Z_n^c} \pi_*. $$
It is easy to verify that $\pi_{*n} \in \Pi(\mu, \nu)$. We have $$ \begin{align} I[\pi_*] &\le I[\pi_{*n}] \\ &= \int_{X\times Y} c d \pi_{*n} \\ &= \pi_*(Z_n)\int_{Z_n} c d \tilde \pi_n + \int_{Z_n^c} c d \pi_* \\ &\le \pi_*(Z_n) I_n [\tilde \pi_n] + 2 \delta \|c\|_\infty \\ & \le I_n [\tilde \pi_n] + 2\delta \|c\|_\infty. \end{align} $$
Now we introduce $$ J_n: L^{1}(\mu_{n}) \times L^{1}(\nu_{n}) \to \mathbb R, (\varphi, \psi) \mapsto \int_{X_{n}} \varphi d \mu_{n}+\int_{Y_{n}} \psi d \nu_{n}. $$
Let $\Phi_{c,n}$ be the set of all $(\varphi, \psi) \in L^{1}(\mu_n) \times L^{1}(\nu_n)$ satisfying $\varphi(x)+\psi(y) \leq c(x, y)$ for $\mu_n$-almost all $x \in X_n$ and $\nu_n$-almost all $y \in Y_n$. By step 1, we know that $$ \inf_{\pi \in \Pi(\mu_n, \nu_n)} I_n [\pi] = \sup_{(\varphi, \psi) \in \Phi_{c,n}} J_n (\varphi, \psi). $$
In particular, there is $(\varphi_n, \psi_n) \in \Phi_{c,n}$ such that $$ J_n (\varphi_n, \psi_n) \ge \sup_{(\varphi, \psi) \in \Phi_{c,n}} J_n (\varphi, \psi) - \delta. $$
To ensure $\varphi_n (x)+\psi_n (y) \leq c(x, y)$ for all $x \in X_n$ and all $y \in Y_n$. We redefine $\varphi_n, \psi_n$ such that $\varphi_n, \psi_n$ take value $-\infty$ on their null sets respectively. WLOG, we assume $\delta <1$. Then $J_n (\varphi_n, \psi_n) \ge -1$. Then there is $(x_n, y_n) \in Z_n$ such that $$ \varphi_n(x_n)+\psi_n(y_n) \ge -1. $$
If we replace $(\varphi_n, \psi_n)$ by $(\varphi_n +s, \psi_n-s)$ for some real number $s$, we do not change the value of $J_n (\varphi_n, \psi_n)$, so the resulting couple is still admissible. By a proper choice of $s$, we can ensure $$ \varphi_n (x_n) \ge -\frac{1}{2} \quad \text{and} \quad\psi_n (y_n) \ge -\frac{1}{2}. $$
This implies for all $(x,y) \in Z_n$, $$ \varphi_n(x) \le c(x, y_n) - \psi_n(y_n) \le c(x, y_n) + 1/2,\\ \psi_n(y) \le c(x_n, y) - \varphi_n(x_n) \le c(x_n, y) + 1/2. $$
Define $\psi_n^c$ by $$ \psi_n^c (x) := \inf_{y\in Y_n} [c(x, y) - \psi_n(y)] \quad \forall x \in X_n. $$
Moreover, $$ \inf_{y\in Y_n} [c(x,y)-c(x_n, y)] - 1/2 \le \psi_n^c (x) \le c (x, y_n) - \psi_n (y_n) \le c(x, y_n) + 1/2. $$
It follows from $c$ is bounded that $\psi_n^c$ is bounded. Hence $\psi_n^c \in L^1(\mu)$. It follows from $\varphi_n(x)+\psi_n(y) \leq c(x, y)$ that $\varphi_n \le \psi_n^c$ on $X_n$. This implies $J_n (\psi_n^c , \psi_n) \ge J_n (\varphi_n, \psi_n)$. Similarly, we define $\psi_n^{cc}$ by $$ \psi_n^{cc} (y) := \inf_{x\in X_n} [c(x, y) - \psi_n^c (x)] \quad \forall y \in Y_n. $$
For all $y\in Y_n$, we have $$ \inf_{x\in X_n} [c(x,y)-c(x, y_n)] - 1/2 \le \psi_n^{cc} (y) \le c(x_n, y) - \psi_n^c (x_n) \le c(x_n, y) - \varphi_n (x_n) \le c(x_n,y)+1/2. $$
So $\psi_n^{cc}$ is also bounded and thus $\psi_n^{cc} \in L^1(\nu)$. Moreover, $$ \psi_n^c (x) + \psi_n^{cc} (y) = \psi_n^c (x) + \inf_{z \in X} [c(z, y) - \psi_n^c(z)] \le \psi_n^c(x) + [c(x, y) - \psi_n^c(x)] = c(x,y). $$
Then $(\psi_n^c, \psi_n^{cc}) \in \Phi_c$. Also $$ \psi_n^{cc} (y) = \inf_{x\in X_n} [c(x, y) - \psi_n^c(x)] \ge \inf_{x\in X_n} [c(x, y) - [c(c,y) - \psi_n (y)]] = \psi_n (y) \quad \forall y \in Y_n. $$
This implies $$ J_n (\psi_n^c , \psi_n^{cc}) \ge J_n(\psi_n^c, \psi_n) \ge J_n (\varphi_n, \psi_n). $$
In particular, $$ \psi_n^c (x) \ge - \|c\|_\infty -1/2 \quad \text{and} \quad \psi_n^{cc} (y) \ge - \|c\|_\infty -1/2. $$
Indeed, $$ \begin{align*} J(\psi_n^c, \psi_n^{cc}) &= \int_{X} \psi_n^c d \mu+\int_{Y} \psi_n^{cc} d \nu \\ &= \int_{X \times Y}[\psi_n^c (x) + \psi_n^{cc} (y)] d \pi_{*n}(x, y) \\ &= \pi_{*}[Z_n] \int_{Z_n}[ \psi_n^c (x) + \psi_n^{cc} (y)] d \tilde \pi_{n}(x, y) +\int_{Z_n^{c}}[\psi_n^c(x)+\psi_n^{cc}(y)] d \pi_{*}(x, y) \\ &\ge (1-2 \delta) \left (\int_{X_n} \psi_n^c d \mu_n+\int_{Y_n} \psi_n^{cc} d \nu_n \right )-(2\|c\|_{\infty}+1) \pi_{*}[Z_n^c] \\ &\geq(1-2 \delta) J_n( \psi_n^c , \psi_n^{cc} )-2(2\|c\|_{\infty}+1) \delta \\ &\geq(1-2 \delta) J_n({\varphi}_n, {\psi}_n)-2(2\|c\|_{\infty}+1) \delta \\ &\geq(1-2 \delta) \left ( \inf_{\pi \in \Pi(\mu_n, \nu_n)} I_n [\pi] - \delta \right)-2(2\|c\|_{\infty}+1) \delta \\ &= (1-2 \delta)( I_n [\tilde \pi_n] -\delta)-2(2\|c\|_{\infty}+1) \delta \\ &\ge (1-2 \delta)(I[\pi_{*}] -(2\|c\|_{\infty}+1) \delta)-2(2\|c\|_{\infty}+1) \delta. \end{align*} $$
The result then follows by taking the limit $\delta \to 0^+$.