Is maximum likelihood estimation practical through backpropagation?

Question

I am trying to implement MLE using pytorch. I created a simple linear regression data set $\mathcal{D} :=\{X, Y\} = \{(x_1, y_1), ...,(x_N, y_N)\}$, so one can write the likelihood function:

$$ p(y|x) = \mathcal{N}(y|f(\boldsymbol x),\,\sigma^{2}) \ \text{where} \, \boldsymbol x \in \mathbb{R}^D ,y \in \mathbb{R}$$ $$f(x) = \boldsymbol x^T \boldsymbol \theta $$

Assuming that all observations are $\text{i.i.d.}$, likelihood factorizes to: $$ p(Y|X, \theta)=p(y_1,...,y_N|\boldsymbol x_1, ... , \boldsymbol x_N, \boldsymbol \theta) = \prod_{n=1}^{N} p(y_n|\boldsymbol x_n, \boldsymbol \theta) = \prod_{n=1}^{N} \mathcal{N}(y_n|f(\boldsymbol x_n),\,\sigma^{2})$$

And our objective is: $$ \boldsymbol \theta_{ML} = \text{argmax} \ p(Y|X, \boldsymbol \theta)$$

Now, here comes my problem, when I initialize $\boldsymbol \theta$ and calculate the likelihood by multiplying the density function of each observation I end up with a very small number that eventually becomes zero (especially if $N$ is large), and I can't even call tensor.backward() to use gradient descent to minimize the minus log of the likelihood function.

However, for some cases I managed to get a working solution by multiplying the density function of each observation by a constant to avoid ending up with zero, but it doesn't always work.

Note: I'm aware of the closed-form solution for the optimum $\boldsymbol \theta^*$, I am just trying to reach optimum parameters through gradient descent and backpropagation of negative log likelihood.

So, is there a work around for calculating the likelihood of the training data set without ending up with zero?

Answer 1

Traditionally, machine learning approaches maximize the log-likelihood instead of the actual likelihood. That is, instead of maximizing $p(Y|X, \theta)$ wrt $\theta$, we instead maximize

$$\log p(Y | X, \theta) = \log \prod_i p(y_i | x_i, \theta) = \sum_i \log p(y_i | x_i, \theta).$$

Because the logarithm turns products into summation, the computation of the log likelihood is a sum of large negative numbers, as opposed to a product of small positive numbers. This leads to better numerical stability compared to the actual likelihood, as you have noticed. Equivalently, you can also minimize the negative of the log likelihood, which is exactly the same as maximizing the log likelihood.

From your post, I'm guessing that what you are trying to do is roughly (in pseudocode):

likelihood = 1
for datapoint in dataset:
    likelihood *= density(datapoint)
log_likelihood = log(likelihood)
log_likelihood.backward()

what you should be doing instead is:

log_likelihood = 0
for datapoint in dataset:
    log_likelihood += log(density(datapoint))
log_likelihood.backward()