Let me address various points of the discussions regarding the use of scoring rules. I think the argument for using (strictly) proper scoring rules in multi-label classification is the same as in binary or multi-class classification. They allow us to compare the quality of statistical information on the labeling provided by forecasters or models.
But how do strictly proper scoring rules apply to multi-label classification?
There are various ways of application, and to explain the two which I find most insightful I borrow the two major "problem transformation methods" for multi-label classification which are mentioned in the corresponding Wikipedia article.
1) Transformation into multi-class classification
This essentially means that we do classification on the label powerset, i.e. instead of looking at $k$ labels we consider $2^k$ categories which reflect all possible label combinations. For instance, if the labels can be ${a,b}$ we would use the categories $\{\emptyset, \{a\}, \{b\}, \{a,b\} \}$. We could then require multi-label classifications to come in the form of probability distributions on the label powerset. See also Ben's answer to What is the statistical model for a multi-label problem?. These probability distributions can be evaluated via proper scoring rules:
Let $m=2^k$ be the size of the label powerset and assume $p = (p_1, \ldots, p_m)$ are probabilities for the different label combinations. Let $l$ be a labelling and $i(l)$ the index of the label combination which corresponds to $l$ (so in our example above $i(\{a\}) = 2$). Two popular examples:
The multi-label log score / cross entropy loss is
$$
S_{\log} (p, l) = - \log ( p_{i(l)} )
$$
If we have $n$ datapoints $l_1, \ldots, l_n$, then using the average $S_{\log}$ of the data can be interpreted as the log-likelihood of a multinomial distribution with $n$ trials, $m=2^k$ possible events and probabilities $p_1, \ldots, p_m$.
Alternatively, the multi-label probability/quadratic score is
$$
S_q (p, l) = \sum_{j=1}^{m} ( p_j - \mathbb{I} \{j = i(l) \} )^2
$$
Both are strictly proper scoring rules, since they are just known scoring rules applied to a certain discrete distribution.
2) Transformation into binary classifications
In this case we have simply $k$ probability forecasting problems which share one observation. Hence, we can apply (strictly) proper scoring rules to each individual label forecast. Assume $p= (p_1, \ldots, p_k)$, where $p_i$ is the probability that label $i$ is present (The components of $p$ do not sum to 1!) Let $S: [0,1] \times {0,1} \to \mathbb{R}$ be a (strictly) proper scoring rule for probabilities and let $y(l)$ be the vector which indicates which label is present (so in our example above $y(\{a\}) = (1, 0)$ and $y(\{a,b\}) = (1,1)$.) Then
$$
S_\mathrm{bin} (p, l) = \sum_{j=1}^{k} S(p_j, (y(l))_j )
$$
is a proper scoring rule. For $S$ we could for example use the log score or the Brier score.
Using this evaluation approach does not imply that we assume the labels occur independently. It means that we ignore the information on label dependence in the evaluation. Thinking about it in the sense of full probability distributions as discussed in 1) we only evaluate the marginals of the full distribution. As a result, a model which gets the marginals right, but dependence wrong, can achieve the same score as the best possible model which takes independence into account.
Ignoring dependence in the evaluation ignores information, but it has the advantage that it simplifies the evaluation. In particular, it allows for easy assessment of calibration and computation of score decompositions, since these are well-known for binary problems.
In comparison to a standard multi-class setting, I see no drawbacks or advantages of proper scoring rules. The only difference is the size of the possible observations ($k$ vs $2^k$), which can make the evaluation expensive to compute.
