When dealing with non-binary discrete-choice outcomes, one common way of modeling such problems as a function of some covariates is through a multinomial logit/logistic model, in which there is one set of coefficients per class/choice, each producing a linear score for it, and the probability of a given observation belonging to each class is calculated as: $$ P(y = k | \mathbf{x}) = \frac{ \text{exp}( \mathbf{x} \mathbf{\beta}_k ) }{ \sum_{i=1}^n \text{exp}( \mathbf{x} \mathbf{\beta}_i ) } $$
Given this definition of $P(y = k)$, it would be my intuition that this model should be fit by maximizing the log-likelihood of the observed outcome for each row - i.e.: $$ \sum_{i=1}^m \text{log}( P(y=y_m | \mathbf{x}_m) ) $$
However, if I look at popular software for fitting such models such as GLMNET or scikit-learn, I see that they instead calculate their optimization target as a sum over all the possible values of $y$ - i.e.: $$ \sum_{i=1}^m ( \sum_{k=1}^K I(y_i = k) \times \mathbf{x} \mathbf{\beta}_k - \text{log}( \sum_{j} \text{exp}( \mathbf{x} \mathbf{\beta}_j ) ) ) $$
(where $I(.)$ is the indicator function taking a value of $1$ when the argument is true and zero otherwise)
My question is: why would calculation of the likelihood need to include all the possible values of the response variable if the probability for only the observed outcome of each row already incorporates the coefficients from all the choices?
References:
https://glmnet.stanford.edu/articles/glmnet.html#multinomial-regression-family-multinomial-