  Kernel techniques are among the most popular and powerful approaches of data
science. Among the key features that make kernels ubiquitous are (i) the number
of domains they have been designed for, (ii) the Hilbert structure of the
function class associated to kernels facilitating their statistical analysis,
and (iii) their ability to represent probability distributions without loss of
information. These properties give rise to the immense success of
Hilbert-Schmidt independence criterion (HSIC) which is able to capture joint
independence of random variables under mild conditions, and permits closed-form
estimators with quadratic computational complexity (w.r.t. the sample size). In
order to alleviate the quadratic computational bottleneck in large-scale
applications, multiple HSIC approximations have been proposed, however these
estimators are restricted to $M=2$ random variables, do not extend naturally to
the $M\ge 2$ case, and lack theoretical guarantees. In this work, we propose an
alternative Nystr\"om-based HSIC estimator which handles the $M\ge 2$ case,
prove its consistency, and demonstrate its applicability in multiple contexts,
including synthetic examples, dependency testing of media annotations, and
causal discovery.
