Are there any general strategies to prove $K(x,y)$ is a machine learning kernel? (I.e. always defines a covariance matrix)?

So there are certain functions of two variables such as the standard Gaussian/radial function $K(x_i,x_j) = e^{-(x_i-x_j)^2}$ which are "kernels" as machine learning calls them, meaning that for any finite set of $N$ points $\{x_i\}$, the $N \times N$ matrix whose $(i,j)$th entry is $K(x_i,x_j)$ is a covariance matrix, meaning symmetric positive-semidefinite. Here $N$ can be any positive integer, and the $x_i$ can be any valid inputs (e.g. any real numbers in the case of the standard Gaussian/radial function).

So, given a function $K(x,y)$, are there any general strategies (or general function properties, such as symmetric and concave/convex) so that it be verified that $K$ is a kernel when there are infinitely many (in fact, even uncountably) matrices that can be formed from $K$ which all need to be covariance matrices?

First, I'll list some characterizations, and then make some remarks on how all this can actually be used.

The definition
Let $X$ be some set. We call a function $k: X\times X \to \mathbb{R}$ a kernel, if there is some Hilbert space $(H,<\bullet,\bullet>)$ and a function $\phi: X \to H$ such that $k(x,y) = <\phi(x),\phi(y)>$ holds for all $x,y\in X$.

For the most basic cases this definition can be verified directly. A trivial fact that is useful to get started is that the ordinary scalar product on $\mathbb{R}$ is a kernel, and that for example $(x,y) \mapsto x_i y_i$ for a fixed index $i$ is a kernel.

Characterization with Gram's matrices
A function $k: X \times X \to \mathbb{R}$ is a kernel iff for each choice of points $\{x_1, ..., x_N\}$ from $X$ the matrix $K = (k(x_i, x_j))_{i,j}$ is positive semidefinite.

This characterization becomes very useful if we want to prove that certain compositions of kernels are again kernels, because all we have to do is to verify a simple property of a finite matrix. In particular, this criterion is used in the derivation of the following rules:

Composing kernels
Suppose $k, m$ are kernels. Then the following functions are also kernels:

1. $k + m$
2. $c k$ for positive constants $c$
3. $k \cdot m$ (pointwise product)
4. $p\circ k$ for arbitrary polynomials with positive coefficients $p$, for example $(x,y) \mapsto 777 \cdot k(x,y)^7 + 42\cdot k(x,y) ^{123}$.

Now the kernel-property also carries over on certain functions that can be approximated by polynomials, for example:

1. $(x,y) \mapsto exp(k(x,y))$ is again a kernel.

Mercer's condition
If $X$ is equipped with a measure $\mu$ it's enough to verify that for all square integrable $f\in L^2(X,\mu)$ it holds:

$\qquad \int \int f(x)k(x,y)f(y) \mu[dx] \mu[dy] \geq 0$

which is particularly useful for all kind of discrete counting measures.

What usually happens
Usually one starts with the definition, then one proves the characterization with Gram's matrices, then uses the Gram matrices in order to prove that we can compose kernels in certain ways, and from then on one more or less forgets all the scary details, and starts to build kernels using the rules of composition.

The situation is similar to what happens with continuous functions: one shows that some basic functions (e.g. identity function on $\mathbb{R}$) is continuous, then one shows that sums and products and compositions are continuous, and that the continuity carries over to certain limit functions. From then on, one never uses any $\epsilon-\delta$ arguments to show that $37 x ^{-17}+\arctan\left(\frac{\sqrt{x^137+\Gamma(59\cdot x)}}{log(1+x^2)}\right)$ is continuous. It the same with kernels: one shows a few basic rules for building more complex kernels from simpler ones, and then uses these basic rules to compose complex problem-specific kernels.