Given a function \(f\) on the positive half-line \(\R_+\) and a sequence (finite or infinite) of points \(X=\{x_k\}_{k=1}^\omega\) in \(\R^n\), we define and study matrices \(\kS_X(f)=\|f(|x_i-x_j|)\|_{i,j=1}^\omega\) called Schoenberg's matrices. We are primarily interested in those matrices which generate bounded and invertible linear operators \(S_X(f)\) on \(\ell^2(\N)\). We provide conditions on \(X\) and \(f\) for the latter to hold. If \(f\) is an \(\ell^2\)-positive definite function, such conditions are given in terms of the Schoenberg measure \(\sigma(f)\). We also approach Schoenberg's matrices from the viewpoint of harmonic analysis on \(\R^n\), wherein the notion of the strong \(X\)-positive definiteness plays a key role. In particular, we prove that \emph{each radial \(\ell^2\)-positive definite function is strongly \(X\)-positive definite} whenever \(X\) is separated. We also implement a "grammization" procedure for certain positive definite Schoenberg's matrices. This leads to Riesz--Fischer and Riesz sequences (Riesz bases in their linear span) of the form \(\kF_X(f)=\{f(x-x_j)\}_{x_j\in X}\) for certain radial functions \(f\in L^2(\R^n)\). Examples of Schoenberg's operators with various spectral properties are presented.