  Kernel quadratures and other kernel-based approximation methods typically
suffer from prohibitive cubic time and quadratic space complexity in the number
of function evaluations. The problem arises because a system of linear
equations needs to be solved. In this article we show that the weights of a
kernel quadrature rule can be computed efficiently and exactly for up to tens
of millions of nodes if the kernel, integration domain, and measure are fully
symmetric and the node set is a union of fully symmetric sets. This is based on
the observations that in such a setting there are only as many distinct weights
as there are fully symmetric sets and that these weights can be solved from a
linear system of equations constructed out of row sums of certain submatrices
of the full kernel matrix. We present several numerical examples that show
feasibility, both for a large number of nodes and in high dimensions, of the
developed fully symmetric kernel quadrature rules. Most prominent of the fully
symmetric kernel quadrature rules we propose are those that use sparse grids.
