How to set up neural network to output ordinal data? I have a neural network set up to predict something where the output variable is ordinal.  I will describe below using three possible outputs A < B < C.  
It is pretty obvious how to use a neural network to output categorical data: the output is just a softmax of the last (usually fully connected) layer, one per category, and the predicted category is the one with the largest output value (this is the default in many popular models).  I have been using the same setup for ordinal values.  However, in this case the outputs often don't make sense, for example the network outputs for A and C are high but B is low: this is not plausible for ordinal values.
I have one idea for this, which is to calculate loss based on comparing the outputs with 1 0 0 for A, 1 1 0 for B, and 1 1 1 for C.  The exact thresholds can be tuned later using another classifier (eg Bayesian) but this seems to capture the essential idea of an ordering of inputs, without prescribing any specific interval scale.
What is the standard way of solving this problem?  Is there any research or references that describe the pros and cons of different approaches?
 A: Update: Meanwhile, I delved more into the topic and even wrote a package implementing many ordinal losses from the literature. It includes the loss I mention here (ordinal encoding) but many others as well.
I believe what most people do is to simply treat ordinal classification as a generic multi-class classification. So, if they have $K$ classes, they will have $K$ outputs, and simply use cross-entropy as the loss.
But some people have managed to invent a clever encoding for your ordinal classes (see this stackoverflow answer). It's a sort of one-hot encoding,

*

*class 1 is represented as [0 0 0 0 ...]


*class 2 is represented as [1 0 0 0 ...]


*class 3 is represented as [1 1 0 0 ...]
i.e. each neuron is predicting the probability $P(\hat y < k)$. You still have to use a sigmoid as the activation function, but I think this helps the network understanding some continuity between classes, I don't know. Afterwards, you do a post-processing (np.sum) to convert the binary output into your classes.
This strategy resembles the ensemble from Frank and Hall, and I think this is the first publication of such.
A: I think the approach to only encode the ordinal labels as

*

*class 1 is represented as [0 0 0 0 ...]


*class 2 is represented as [1 0 0 0 ...]


*class 3 is represented as [1 1 0 0 ...]
and use binary cross-entropy as the loss function is suboptimal. As mentioned in the comments, it might happen that the predicted vector is for example [1 0 1 0 ...]. This is undesirable for making predictions.
The paper Rank-consistent ordinal regression for neural networks describes how to restrict the neural network to make rank-consistent predictions. You have to make sure that the last layer shares its weights, but should have different biases. You can implement this in Tensorflow by adding the following as the last part of the network (credits for https://stackoverflow.com/questions/59656313/how-to-share-weights-and-not-biases-in-keras-dense-layers):
class BiasLayer(tf.keras.layers.Layer):
    def __init__(self, units, *args, **kwargs):
        super(BiasLayer, self).__init__(*args, **kwargs)
        self.bias = self.add_weight('bias',
                                    shape=[units],
                                    initializer='zeros',
                                    trainable=True)

    def call(self, x):
        return x + self.bias


# Add the following as the output of the Sequential model
model.add(keras.layers.Dense(1, use_bias=False))
model.add(BiasLayer(4))
model.add(keras.layers.Activation("sigmoid"))

Note that the number of ordinal classes here is 5, hence the $K-1$ biases.
I tested the difference in performance on actual data, and the predictive accuracy improved substantially. Hope this helps.
