Following the line of thought by P. Bouboulis & S. Theodoridis, [6], we take up a program of recovering kernel methods (as employed in signal analysis and machine learning theory) from real RKHS and kernels, to the complex domain. We solve the maximum problem \(\sup \, \big \{ \sum _{j=1}^p \big | f( z_j ) \big |^2 \,: \, \Vert f \Vert ^2 \le E \big \}\) in the complex RKHS of holomorphic \(L^2\) functions \(f: \Omega \rightarrow {\mathbb {C}}\) , for any bounded domain \(\Omega \subset {\mathbb {C}}^n\) and any finite set of points \(z_1 \,, \, \cdots \,, \, z_p \in \Omega \) , and apply the result to the space \(L^2 H({\mathbb {B}}^n )\) of holomorphic \(L^2\) functions on the unit ball \({\mathbb {B}}^n \subset {\mathbb {C}}^n\) . We produce sampling expansions of functions \(f \in L^2 H(\Omega )\) associated to infinite sequences \(\{ \zeta _k \}_{k \ge 0} \subset \Omega \) , by starting from complete orthonormal systems \(\{ \phi _\nu \}_{\nu \ge 0} \subset L^2 H(\Omega )\) and approximating each \(\phi _\nu \) uniformly on \(\Omega \) by a linear combination of reproducing kernels. The means to said approximation are provided by the Faber-Kaczmarz-Mycielski algorithm \({\mathscr {A}}(h)\) learning (cf. [23]) from the data \(\big \{ \big ( \zeta _k \,, \, \phi _\nu (\zeta _k ) \big ) \big \}_{k \ge 0}\) and producing an approximating sequence \(\big \{ {\mathscr {J}}_k \, \phi _\nu \big \}_{k \ge 0} \subset L^2 H(\Omega )\) .