Solution (source code)

= Solution

Let $B$ be the integral operator with kernel $\beta$. Since $Y=BX+\varepsilon$ and the centered error is independent of $X$,
$$
\mathbb E[a_kb_l]
=\mathbb E\!\left[a_k\langle BX,u_l\rangle\right].
$$
Writing $\beta_{lk}=\langle B\phi_k,u_l\rangle$ and using $\mathbb E[a_ka_j]=\lambda_k\mathbf1_{\{j=k\}}$ gives
$$
\mathbb E[a_kb_l]=\lambda_k\beta_{lk}.
$$
Expanding the <function-on-function linear model> kernel in the product basis therefore yields
$$
\boxed{\beta(t,s)=
\sum_{k=1}^\infty\sum_{l=1}^\infty
\frac{\mathbb E[a_kb_l]}{\lambda_k}\phi_k(s)u_l(t)}.
$$