## I.1 Closed-form expression for $\hat{h}_{\omega}$

Let $\omega \in \mathbb{R}^d$. By the first-order optimality condition, the gradient of $\hat{L}_{in}$ with respect to its second argument must vanish at $\hat{h}_{\omega}$, i.e., $\partial_h \hat{L}_{in}(\omega, \hat{h}_{\omega}) = 0$. Proposition B.1 implies that $\hat{h}_{\omega}$ satisfies the following equation:

$$ \frac{1}{n} \sum_{i=1}^{n} \left[ (\hat{h}_{\omega}(x_i) - \omega^{\top} \phi(t_i)) K(x_i, \cdot) \right] + \lambda \hat{h}_{\omega} = 0, $$

with $(x_i, t_i)_{1 \le i \le n}$ being $n$ samples drawn from the distribution $\mathbb{P}$. After using the reproducing property of the RKHS $\mathcal{H}$ and rearranging the terms, we arrive at the following closed-form expression for $\hat{h}_{\omega}$:

$$ \hat{h}_{\omega} = (\hat{\Sigma}_{\lambda}^{-1} \hat{\Phi})^{\top} \omega \in \mathcal{H}, $$

where $\hat{\Sigma}_{\lambda} = \hat{\Sigma} + \lambda \text{Id}_{\mathcal{H}}$ is an operator from $\mathcal{H}$ to $\mathcal{H}$ with $\hat{\Sigma} = \frac{1}{n} \sum_{i=1}^{n} K(x_i, \cdot) \otimes K(x_i, \cdot)$ being the empirical covariance operator and $\hat{\Phi} = \frac{1}{n} \sum_{i=1}^{n} \phi(t_i) K(x_i, \cdot) = (\hat{\Phi}_1, \dots, \hat{\Phi}_d)^{\top} \in \mathcal{H}^d$. Next, we compute a closed-form expression for $\hat{\Sigma}_{\lambda}^{-1} \hat{\Phi}$, which can be determined as the solution $\hat{b} = (\hat{b}_1, \dots, \hat{b}_d)^{\top} \in \mathcal{H}^d$ of the following minimization problem:

$$ \hat{b}_l = \arg \min_{b_l \in \mathcal{H}} \frac{1}{2} b_l^{\top} \hat{\Sigma}_{\lambda} b_l - b_l^{\top} \hat{\Phi}_l, \quad \text{for any } 1 \le l \le d. $$

After expanding the terms and rearranging, this minimization problem is equivalent to:

$$ \hat{b}_l = \arg \min_{b_l \in \mathcal{H}} \Psi_l(b_l(x_1), \dots, b_l(x_n), \|b_l\|_{\mathcal{H}}), \quad \text{for any } 1 \le l \le d, $$

where, for any $e_1, \dots, e_n, e \in \mathbb{R}$, $\Psi_l(e_1, \dots, e_n, e) = \frac{1}{2n} \sum_{i=1}^{n} e_i^2 - \frac{1}{n} \sum_{i=1}^{n} \phi_l(t_i) e_i + \frac{\lambda}{2} e^2$. By the representer theorem, for any $1 \le l \le d$, $\hat{b}_l$ can be expressed as:

$$ \hat{b}_l = \sum_{i=1}^{n} \hat{c}_{i,l} K(x_i, \cdot), $$

where $\hat{c}_l = (\hat{c}_{1,l}, \dots, \hat{c}_{n,l})^{\top} \in \mathbb{R}^n$ satisfies:

$$ \hat{c}_l = \arg \min_{\mathbf{c}_l \in \mathbb{R}^n} \Psi_l ([\mathbf{K} \mathbf{c}_l]_1, \dots, [\mathbf{K} \mathbf{c}_l]_n, \mathbf{c}_l^{\top} \mathbf{K} \mathbf{c}_l) := \frac{1}{2n} \mathbf{c}_l^{\top} \mathbf{K}^2 \mathbf{c}_l - \frac{1}{n} \mathbf{F}_l^{\top} \mathbf{K} \mathbf{c}_l + \frac{\lambda}{2} \mathbf{c}_l^{\top} \mathbf{K} \mathbf{c}_l, $$

where $\mathbf{F}_l = (\phi_l(t_1), \dots, \phi_l(t_n))^{\top} \in \mathbb{R}^n$. By the first-order optimality condition, we have:

$$ \nabla_{\mathbf{c}_l} \Psi_l ([\mathbf{K} \hat{c}_l]_1, \dots, [\mathbf{K} \hat{c}_l]_n, \hat{c}_l^{\top} \mathbf{K} \hat{c}_l) = 0, \text{ which results in } \hat{c}_l = (\mathbf{K} + n\lambda \mathbf{1}_{n \times n})^{-1} \mathbf{F}_l \in \mathbb{R}^n. $$

Using this, we obtain:

$$ \hat{b}_l = \hat{c}_l^{\top} (K(x_1, \cdot), \dots, K(x_n, \cdot))^{\top}, \text{ for any } 1 \le l \le d, \text{ and thus: } \hat{h}_{\omega} = \hat{b}^{\top} \omega. $$

Now that we have obtained a closed-form expression for $\hat{h}_{\omega}$, we can express $\hat{\mathcal{F}}(\omega)$ and $\nabla \mathcal{F}(\omega)$ in closed-form, as we will see next.

## I.2 Plug-in estimators for $\mathcal{F}(\omega)$ and $\nabla \mathcal{F}(\omega)$

Let $(\tilde{x}_j, \tilde{y}_j)_{1 \le j \le m}$ be $m$ samples drawn from $\mathbb{Q}$ and $\omega \in \mathbb{R}^d$. We have:

$$ \hat{\mathcal{F}}(\omega) = \frac{1}{2m} \sum_{j=1}^{m} (\hat{h}_{\omega}(\tilde{x}_j) - \tilde{y}_j)^2 = \frac{1}{2m} \|\hat{\mathbf{B}}\omega - \tilde{\mathbf{y}}\|^2, $$

where $\hat{\mathbf{B}} = [\hat{b}(\tilde{x}_1), \dots, \hat{b}(\tilde{x}_m)]^{\top} \in \mathbb{R}^{m \times d}$ and $\tilde{\mathbf{y}} = (\tilde{y}_1, \dots, \tilde{y}_m)^{\top} \in \mathbb{R}^m$. For any $1 \le l \le d$ and $1 \le j \le m$, we have:

$$ \hat{b}_l(\tilde{x}_j) = [\mathbf{F}_l]^{\top} (\mathbf{K} + n\lambda \mathbf{1}_{n \times n})^{-1} [\bar{\mathbf{K}}^{\top}]_j. $$
As a consequence, we obtain:

$$ \hat{b}(\tilde{x}_j) = (\hat{b}_1(\tilde{x}_j), \dots, \hat{b}_d(\tilde{x}_j))^{\top} = \hat{\mathbf{C}}^{\top} [\overline{\mathbf{K}}^{\top}]_j \in \mathbb{R}^d, $$

where $\hat{\mathbf{C}} = [\hat{c}_1, \dots, \hat{c}_d] = (\mathbf{K} + n\lambda \mathbf{1}_{n \times n})^{-1} \mathbf{F} \in \mathbb{R}^{n \times d}$, with $\mathbf{F} = [\mathbf{F}_1, \dots, \mathbf{F}_d] \in \mathbb{R}^{n \times d}$. This implies that $\hat{\mathbf{B}} = \overline{\mathbf{K}} \hat{\mathbf{C}} \in \mathbb{R}^{m \times d}$, and hence:

$$ \hat{\mathcal{F}}(\omega) = \frac{1}{2m} \|\overline{\mathbf{K}} \hat{\mathbf{C}}\omega - \tilde{\mathbf{y}}\|^2 = \frac{1}{2m} \omega^{\top} (\overline{\mathbf{K}} \hat{\mathbf{C}})^{\top} \overline{\mathbf{K}} \hat{\mathbf{C}} \omega - \frac{1}{m} \tilde{\mathbf{y}}^{\top} \overline{\mathbf{K}} \hat{\mathbf{C}} \omega + \frac{1}{2m} \|\tilde{\mathbf{y}}\|^2. \quad (29) $$

On the other hand, using Appendix C, we get:

$$ \overline{\nabla \mathcal{F}}(\omega) = \frac{1}{m} \hat{\mathbf{C}}^{\top} \overline{\mathbf{K}}^{\top} (\overline{\mathbf{K}} \hat{\mathbf{C}} \omega - \tilde{\mathbf{y}}) = \frac{1}{m} [(\overline{\mathbf{K}} \hat{\mathbf{C}})^{\top} \overline{\mathbf{K}} \hat{\mathbf{C}} \omega - (\overline{\mathbf{K}} \hat{\mathbf{C}})^{\top} \tilde{\mathbf{y}}] \in \mathbb{R}^d. \quad (30) $$

The exact expressions of $\mathcal{F}(\omega)$ and $\overline{\nabla \mathcal{F}}(\omega)$ involve expectations and are therefore intractable to compute analytically. A natural approach is to approximate these quantities using their plug-in estimators $\hat{\mathcal{F}}(\omega)$ and $\overline{\nabla \mathcal{F}}(\omega)$, evaluated with a very large number of inner and outer samples, $n$ and $m$. However, this approach quickly becomes computationally and memory-intensive. In particular, storing the kernel matrices $\mathbf{K}$ and $\overline{\mathbf{K}}$ requires $O(n^2)$ and $O(nm)$ space, respectively. Moreover, computing the inverse $(\mathbf{K} + n\lambda \mathbf{1}_{n \times n})^{-1}$ incurs a cubic time complexity of $O(n^3)$, which is prohibitive for large-scale applications. To alleviate these computational bottlenecks, potential strategies rely on classical techniques in kernel methods such as Random Fourier Features (RFF), which approximate kernel functions in a finite-dimensional feature space and enable more efficient gradient computations [60], and Nyström approximations, which mitigate the computational burden of full kernel matrices by using a low-rank approximation of the kernel [75]. In our experiments, we leverage the closed-form expressions of the plug-in estimators, and replace the kernel evaluations with their approximations via RFF. This enables us to construct efficient and scalable approximations of $\mathcal{F}(\omega)$ and $\overline{\nabla \mathcal{F}}(\omega)$, while significantly reducing both the memory usage and the computational cost. Our approach will be discussed in the following.

## I.3 Scalable approximations for $\mathcal{F}(\omega)$ and $\overline{\nabla \mathcal{F}}(\omega)$ via random Fourier features

Random Fourier Features (RFF) provide a way to approximate shift-invariant kernels (i.e., kernels that satisfy $K(x, x') = G(x - x')$ for some function $G: \mathcal{X} \to \mathbb{R}$) by mapping the data into a randomized feature space. To avoid the high computational burden of building the full kernel matrix from all pairwise kernel evaluations, RFF use a randomized feature map $\psi: \mathcal{X} \to \mathbb{R}^D$, with $D$ being the number of Fourier features, to approximate the kernel as follows:

$$ K(x, x') \approx \psi(x)^{\top} \psi(x'), \quad \text{for any } x, x' \in \mathcal{X}. $$

Now, we derive the expression of the feature map $\psi$. Let $x, x' \in \mathcal{X}$. By Bochner theorem [16], we have that:

$$ K(x, x') = \frac{1}{(2\pi)^p} \mathbb{E}_{\mathbf{w} \sim \hat{\mathcal{G}}(\mathbf{w})} \left[ e^{i \mathbf{w}^{\top} (x-x')} \right] = \frac{1}{(2\pi)^p} \mathbb{E}_{\mathbf{w} \sim \hat{\mathcal{G}}(\mathbf{w})} \left[ \cos(\mathbf{w}^{\top} (x-x')) \right], \quad (31) $$

where $\hat{\mathcal{G}}$ is the Fourier transform of $G$. For any $b \in \mathbb{R}$, the following product-to-sum identity holds:

$$ 2 \cos(\mathbf{w}^{\top} x + b) \cos(\mathbf{w}^{\top} x' + b) = \cos(2b + \mathbf{w}^{\top}(x + x')) + \cos(\mathbf{w}^{\top}(x - x')). $$

In particular, when $b \sim \mathcal{U}(0, 2\pi)$ (the uniform distribution over $[0, 2\pi]$), we get:

$$ \mathbb{E}_{b \sim \mathcal{U}(0,2\pi)} [2 \cos(\mathbf{w}^{\top} x + b) \cos(\mathbf{w}^{\top} x' + b)] = \mathbb{E}_{b \sim \mathcal{U}(0,2\pi)} [\cos(2b + \mathbf{w}^{\top}(x+x'))] \\ + \cos(\mathbf{w}^{\top}(x-x')). $$

However, we have:

$$ \begin{aligned} \mathbb{E}_{b \sim \mathcal{U}(0,2\pi)} [\cos(2b + \mathbf{w}^{\top}(x+x'))] &= \frac{1}{2\pi} \int_{0}^{2\pi} \cos(2b + \mathbf{w}^{\top}(x+x')) \, db \\ &= \frac{1}{4\pi} [\sin(2b + \mathbf{w}^{\top}(x+x'))]_{b=0}^{b=2\pi} = 0. \end{aligned} $$
Thus:

$$ \mathbb{E}_{b \sim \mathcal{U}(0,2\pi)} [2 \cos(\mathbf{w}^{\top} x + b) \cos(\mathbf{w}^{\top} x' + b)] = \cos(\mathbf{w}^{\top}(x - x')). $$

Substituting this back into Equation (31), we arrive at:

$$ K(x, x') = \mathbb{E}_{\mathbf{w} \sim \frac{1}{(2\pi)^p} \hat{G}(\mathbf{w}), b \sim \mathcal{U}(0, 2\pi)} \left[ \sqrt{2} \cos(\mathbf{w}^{\top} x + b) \sqrt{2} \cos(\mathbf{w}^{\top} x' + b) \right]. $$

Using $D$ samples $\mathbf{w}_1, \dots, \mathbf{w}_D \sim \frac{1}{(2\pi)^p} \hat{G}(\mathbf{w})$ and $b_1, \dots, b_D \sim \mathcal{U}(0, 2\pi)$, we obtain by Monte Carlo estimation:

$$ K(x, x') \approx \sum_{i=1}^{D} \left( \sqrt{\frac{2}{D}} \cos(\mathbf{w}_i^{\top} x + b_i) \right) \left( \sqrt{\frac{2}{D}} \cos(\mathbf{w}_i^{\top} x' + b_i) \right). $$

This implies that:

$$ \psi(x) = \sqrt{\frac{2}{D}} \cos(\mathbf{W}x+b), \text{ where } \mathbf{W} = (\mathbf{w}_1, \dots, \mathbf{w}_D)^{\top} \in \mathbb{R}^{D \times p} \text{ and } b = (b_1, \dots, b_D)^{\top} \in \mathbb{R}^D. $$

In practice, one typically chooses $D \ll n$ and $D \ll m$, which reduces the space complexity of storing $\mathbf{K}$ from $\mathcal{O}(n^2)$ to $\mathcal{O}(nD)$, and that of storing $\bar{\mathbf{K}}$ from $\mathcal{O}(nm)$ to $\mathcal{O}((n+m)D)$. This results in significant computational and memory savings. Using the RFF approach, the two kernel matrices $\mathbf{K}$ and $\bar{\mathbf{K}}$ can then be approximated as:

$$ \mathbf{K} \approx \Xi \Xi^{\top} \text{ and } \bar{\mathbf{K}} \approx \tilde{\Xi} \tilde{\Xi}^{\top}, $$

where $\Xi = [\psi(x_1), \dots, \psi(x_n)]^{\top} \in \mathbb{R}^{n \times D}$ and $\tilde{\Xi} = [\psi(\tilde{x}_1), \dots, \psi(\tilde{x}_m)]^{\top} \in \mathbb{R}^{m \times D}$. A common term in Equations (29) and (30) is $\bar{\mathbf{K}} \hat{\mathbf{C}}$, which can be approximated using the push-through identity as follows:

$$ \bar{\mathbf{K}} \hat{\mathbf{C}} \approx \tilde{\Xi} \Xi^{\top} (\Xi \Xi^{\top} + n \lambda \mathbf{1}_{n \times n})^{-1} \mathbf{F} = \tilde{\Xi} (\Xi^{\top} \Xi + n \lambda \mathbf{1}_{D \times D})^{-1} \Xi^{\top} \mathbf{F} \in \mathbb{R}^{m \times d}. $$

Here, instead of inverting a matrix of size $n \times n$, we invert a matrix of size $D \times D$, which leads to significant computational savings in time, especially when $D \ll n$. Consequently, using this approximation, we get:

$$ \hat{\mathcal{F}}(\omega) \approx \frac{1}{2m} \omega^{\top} \mathbf{J}^{\top} \tilde{\Xi}^{\top} \tilde{\Xi} \mathbf{J} \omega - \frac{1}{m} \omega^{\top} \mathbf{J}^{\top} \tilde{\Xi}^{\top} \tilde{\mathbf{y}} + \frac{1}{2m} \|\tilde{\mathbf{y}}\|^2, $$

$$ \widehat{\nabla \mathcal{F}}(\omega) \approx \frac{1}{m} \left[ \mathbf{J}^{\top} \tilde{\Xi}^{\top} \tilde{\Xi} \mathbf{J} \omega - \mathbf{J}^{\top} \tilde{\Xi}^{\top} \tilde{\mathbf{y}} \right], $$

where $\mathbf{J} = (\Xi^{\top} \Xi + n \lambda \mathbf{1}_{D \times D})^{-1} \Xi^{\top} \mathbf{F} \in \mathbb{R}^{D \times d}$. As mentioned earlier, a very large number of samples $n$ and $m$ is required to obtain accurate approximations of $\mathcal{F}(\omega)$ and $\nabla \mathcal{F}(\omega)$ using the RFF approach. To cope with the issue of storing the two matrices $\Xi \in \mathbb{R}^{n \times D}$ and $\tilde{\Xi} \in \mathbb{R}^{m \times D}$ in memory, we implement this method in blocks. More precisely, we divide our data $(x_i, t_i)_{1 \le i \le n}$ and $(\tilde{x}_j, \tilde{y}_j)_{1 \le j \le m}$ into blocks, then compute $\Xi^{\top} \Xi, \tilde{\Xi}^{\top} \tilde{\Xi}, \Xi^{\top} \mathbf{F}, \tilde{\Xi}^{\top} \tilde{\mathbf{y}}$, and $\|\tilde{\mathbf{y}}\|^2$ as follows:

$$ \Xi^{\top} \Xi = \sum_{i=1}^{n} \psi(x_i) \psi(x_i)^{\top} = \sum_{B \in \mathcal{B}} \sum_{x \in B} \psi(x) \psi(x)^{\top} = \sum_{B \in \mathcal{B}} \Xi_B^{\top} \Xi_B \in \mathbb{R}^{D \times D}, $$

$$ \tilde{\Xi}^{\top} \tilde{\Xi} = \sum_{j=1}^{m} \psi(\tilde{x}_j) \psi(\tilde{x}_j)^{\top} = \sum_{B \in \mathcal{B}} \sum_{\tilde{x} \in B} \psi(\tilde{x}) \psi(\tilde{x})^{\top} = \sum_{B \in \mathcal{B}} \tilde{\Xi}_B^{\top} \tilde{\Xi}_B \in \mathbb{R}^{D \times D}, $$

$$ \Xi^{\top} \mathbf{F} = \sum_{i=1}^{n} \psi(x_i) [\phi_1(t_i), \dots, \phi_d(t_i)] = \sum_{B \in \mathcal{B}} \sum_{(x,t) \in B} \psi(x) [\phi_1(t), \dots, \phi_d(t)] = \sum_{B \in \mathcal{B}} \Xi_B^{\top} \mathbf{F}_B \in \mathbb{R}^{D \times d}, $$

$$ \tilde{\Xi}^{\top} \tilde{\mathbf{y}} = \sum_{j=1}^{m} \psi(\tilde{x}_j) \tilde{y}_j = \sum_{B \in \mathcal{B}} \sum_{(\tilde{x}, \tilde{y}) \in B} \psi(\tilde{x}) \tilde{y} = \sum_{B \in \mathcal{B}} \tilde{\Xi}_B^{\top} \tilde{\mathbf{y}}_B \in \mathbb{R}^D, $$

$$ \|\tilde{\mathbf{y}}\|^2 = \sum_{B \in \mathcal{B}} \|\tilde{\mathbf{y}}\|_B^2, $$

where $\mathcal{B}$ denotes the set of blocks, and the subscript $B$ indicates that the corresponding quantity is computed using only the data contained in block $B$. These block-wise computations make it possible to precisely approximate $\hat{\mathcal{F}}(\omega)$ and $\nabla \hat{\mathcal{F}}(\omega)$ in a scalable manner. As a result, we can efficiently approximate both $\mathcal{F}(\omega)$ and $\nabla \mathcal{F}(\omega)$ through their plug-in estimators when choosing large sample sizes $n$ and $m$.
## I.4 Additional details on the experimental setup

We use the JAX framework [18] to run our experiments on an NVIDIA RTX 6000 ADA GPU. The experiments take approximately 15 hours to complete.

**Choice of the kernel.** In our experiments, we consider the Gaussian kernel defined, for any $x, x' \in \mathcal{X}$, as $K(x, x') = e^{-\frac{\|x-x'\|^2}{2\sigma^2}}$, where $\sigma > 0$ is the bandwidth parameter controlling the smoothness. Since the Gaussian kernel is translation-invariant, Bochner's theorem is applicable. In this case, using the same notations as in Appendix I.3, we have $G(z) = e^{-\frac{\|z\|^2}{2\sigma^2}}$, for any $z \in \mathcal{X}$. Its Fourier transform $\hat{G}$ is then given by $\hat{G}(\mathbf{w}) = (2\pi\sigma^2)^{\frac{p}{2}} e^{-\frac{\sigma^2 \|\mathbf{w}\|^2}{2}}$, for any $\mathbf{w} \in \mathbb{R}^d$. As a consequence, we obtain:

$$ \frac{1}{(2\pi)^p} \hat{G}(\mathbf{w}) = \frac{1}{(2\pi)^p} (2\pi\sigma^2)^{\frac{p}{2}} e^{-\frac{\sigma^2 \|\mathbf{w}\|^2}{2}} = \left(\frac{2\pi}{\sigma^2}\right)^{-\frac{p}{2}} e^{-\frac{\sigma^2 \|\mathbf{w}\|^2}{2}} = \mathcal{N}\left(0, \frac{1}{\sigma^2} \mathbb{1}_{p \times p}\right), $$

which implies that $\mathbf{w}_1, \dots, \mathbf{w}_D \sim \mathcal{N}(0, \frac{1}{\sigma^2} \mathbb{1}_{p \times p})$.

**Choice of the statistical model and hyperparameters.** We set $p = 3, d = 4, \lambda = 0.01$, and $\sigma = 0.2$. We generate synthetic data as follows:

$$ x \sim P_x, \quad t = 2(\mathbb{1}_p^\top x + \epsilon), \quad y = \omega^{\top} \phi(t) + \epsilon, $$

where $\epsilon \sim \mathcal{N}(0, 0.025)$, $\omega^* \sim \mathcal{U}(0, 1)^d$, and $\phi(t) = (\sin(t+1), \dots, \sin(t+d))^\top$. We consider two cases for the distribution $P_x$ of the instrumental variable $x$: (i) a $p$-dimensional standard Gaussian, i.e., $P_x = \mathcal{N}(0, \mathbb{1}_{p \times p})$, and (ii) a $p$-dimensional Student's $t$-distribution with degrees of freedom $\nu \in \{2.1, 2.5, 2.9\}$. All random variables are fixed across runs for reproducibility.

## I.5 Additional experimental results

Here, we retain the same experimental setup as in the main paper and extend the analysis by providing additional experimental results in the scenario where both $m$ and $n$ vary simultaneously over the range 100 to 5000. In Figure 3, we visualize the results using heatmaps for four key quantities: $|\mathcal{F}(\omega_0) - \hat{\mathcal{F}}(\omega_0)|$, $\|\nabla \mathcal{F}(\omega_0) - \nabla \hat{\mathcal{F}}(\omega_0)\|$, $\|\nabla \mathcal{F}(\omega_T)\|$, and $\min_{i=0,...,T} \|\nabla \mathcal{F}(\omega_i)\|$, with $n$ on the $x$-axis and $m$ on the $y$-axis. We report the results only for the case where the instrumental variable $x$ is sampled from a $p$-dimensional standard Gaussian, since the cases where $x$ is sampled from a $p$-dimensional Student's $t$-distribution with degrees of freedom $\nu \in \{2.1, 2.5, 2.9\}$ exhibit similar trends. From the heatmaps, we observe that the lowest errors across all four metrics occur along the diagonal of the plots, i.e., when $m = n$. This pattern suggests that matching the number of samples in the two dimensions leads to more accurate estimation of both the objective function and its gradient, as well as improved convergence behavior during optimization.

<figure><img src="image_5.png" /><figcaption>Figure 3: Illustration of gradient descent on (KBO) for the instrumental variable regression task using synthetic data, with an instrumental variable sampled from a standard Gaussian distribution. The logs of the means of the four quantities across 50 runs are displayed.</figcaption></figure>