@leonardoaraujosantos you could try predicting bit-vectors i.e. create N = np.ceil(log2(num_classes)) = 15 (if num_classes=20k) outputs with each one being a sigmoid. The loss would be per-bit binary cross-entropy and the output probabilities would be converted back to 0/1s using some threshold.
Empirically, this worked well on some experiments with small datasets (MNIST and Fashion-MNIST) but of course, they only have 10 classes.