What's the difference between terms 'link function' and 'canonical link function'? Also, are there any (theoretical) advantages of using one over the other?

For example, a binary response variable can be modeled using many link functions such as logit, probit, etc. But, logit here is considered the "canonical" link function.

The above answers are more intuitive, so I try more rigor.
What is a GLM?
Let $Y=(y,\mathbf{x})$ denote a set of a response $y$ and $p$-dimensional covariate vector $\mathbf{x}=(x_1,\dots,x_p)$ with expected value $E(y)=\mu$. For $i=1,\dots,n$ independent observations, the distribution of each $y_i$ is an exponential family with density
$$
f(y_i;\theta_i,\phi)=\exp\left(\frac{y_i\theta_i-\gamma(\theta_i)}{\phi}+\tau(y_i,\phi)\right) = \alpha(y_i, \phi)\exp\left(\frac{y_i\theta_i-\gamma(\theta_i)}{\phi}\right)
$$
Here, the parameter of interest (natural or canonical parameter) is $\theta_i$, $\phi$ is a scale parameter (known or seen as a nuisance) and $\gamma$ and $\tau$ are known functions. The $n$-dimensional vectors of fixed input values for the $p$ explanatory variables are denoted by $\mathbf{x}_1,\dots,\mathbf{x}_p$. We assume that the input vectors influence (1) only via a linear function, the linear predictor,
$$
\eta_i=\beta_0+\beta_1x_{i1}+\dots+\beta_px_{ip}
$$
upon which $\theta_i$ depends. As it can be shown that $\theta=(\gamma')^{-1}(\mu)$, this dependency is established by connecting the linear predictor $\eta$ and $\theta$ via the mean. More specifically, the mean $\mu$ is seen as an invertible and smooth function of the linear predictor, i.e.
$$
g(\mu)=\eta\ \textrm{or}\ \mu=g^{-1}(\eta)
$$
Now to answer your question:
The function $g(\cdot)$ is called the link function. If the function connects $\mu$, $\eta$ and $\theta$ such that $\eta \equiv\theta$, then this link is called canonical and has the form $g=(\gamma')^{-1}$.
That's it. Then there are a number of desirable statistical properties of using the canonical link, e.g., the sufficient statistic is $X'y$ with
components $\sum_i x_{ij} y_i$ for $j = 1, \dots, p$, the Newton Method and Fisher scoring for finding the ML estimator coincide, these links simplify the derivation of the MLE, they ensure that some properties of linear regression (e.g., the sum of the residuals is 0) hold up or they ensure that $\mu$ stays within the range of the outcome variable.
Hence they tend to be used by default. Note however, that there is no a priori reason why the effects in the model should be additive on the scale given by this or any other link.
