It would be pretty easy to compute vectorized, if only `scipy.stats.rankdata` (and the mstats equivalent, perhaps) supported `axis`.

There is certainly a vectorized solution with `O(n_samples * n_labels ^ 2)` space and time complexity, and there may be something faster but more obfuscated. The former is based on the following:
```python
rank = (y_score[:, :, None] >= y_score[:, None, :]).sum(axis=2)
```
and if you take `(y_score[:, :, None] >= y_score[:, None, :])` and mask it with a broadcasted `y_true`, then its `.sum()` will return `L`. Let me know if that needs more clarification.