Scattered data interpolation aims to reconstruct a continuous (smooth) function that approximates the underlying function by fitting (meshless) data points. There are extensive applications of scattered data interpolation in computer graphics, fluid dynamics, inverse kinematics, machine learning, etc. In this paper, we consider a novel generalized Mercel kernel in the reproducing kernel Banach space for scattered data interpolation. The system of interpolation equations is formulated as a multilinear system with a structural tensor, which is an absolutely and uniformly convergent infinite series of symmetric rank-one tensors. Then we design a fast numerical method for computing the product of the structural tensor and any vector in arbitrary precision. Whereafter, a scalable optimization approach equipped with limited-memory BFGS and Wolfe line-search techniques is customized for solving these multilinear systems. Using the Łojasiewicz inequality, we prove that the proposed scalable optimization approach is a globally convergent algorithm and possesses a linear or sublinear convergence rate. Numerical experiments illustrate that the proposed scalable optimization approach can improve the accuracy of interpolation fitting and computational efficiency.