I am getting the first eigenvalue of a covariance matrix with negative value systematically. Why is this?

by Iago González Rodríguez   Last Updated August 14, 2019 00:19 AM

I am working in a implementation of the Expectation-Maximization algorithm with Missing Data for Mixture of MVNs. You don't have to know anything about this algorithm to help me with my issue.

Let say that my dataset has shape D x N with D = 6.

I compute the estimation of sigma (the covariance matrix). The result is a 6 x 6 matrix (everything ok).

Then I have to compute the PDF of some samples for the distribution with the mean and covariance matrix that I just computed (estimated parameters). In order to compute the PDF, I use scipy.stats.multivariate_normal because it has a pdf() method.

But I always get the following error:

ValueError: the input matrix must be positive semidefinite

This error is due to the covariance matrix. I have checked its eigenvalues, and the first eigenvalue is always a negative value. Some real examples that I got:

import numpy as np

eigvals = np.linalg.eigvals(sigma)  # shape (6,)

Three different results that I got in different executions (the
parameters are randomly generated)

[-406.73080893   94.43623712   57.06170498   73.75668673   69.21443981
[-324.74509361  104.75315794   50.43203113   67.92014983   63.35285505
[-277.14957619   98.17501755   69.8623394    54.49958827   59.4808295

Does anyone know why this is happening?

I am not very familiar with the eigenvalues, so I would greatly appreciate your help.

==================== MORE INFO ====================

Here I include some code of how I compute sigma:

# This is how I compute the covariance matrix
def estimate_sigma_with_nan(xx_est, gamma_k, mu):
    sigma = (xx_est / gamma_k) - np.dot(mu, np.transpose(mu))
    return sigma  # = (d, d)

# This is how I compute xx_est. I will include a image of the maths behind this
exp_prod = np.zeros_like(xx_est_k)
x_i_norm[m] = estimate_x_norm(x_i_norm, mu_k, sigma_k, m, o)
exp_prod[np.ix_(m, m)] = estimate_xx_norm(reshape_(x[:d, i], keep_axes=range(2)), mu_k, sigma_k, m, o)
exp_prod[np.ix_(o, o)] = np.dot(x_i_norm[o], np.transpose(x_i_norm[o]))
exp_prod[np.ix_(o, m)] = np.dot(x_i_norm[o], np.transpose(x_i_norm[m]))
exp_prod[np.ix_(m, o)] = np.dot(x_i_norm[m], np.transpose(x_i_norm[o]))
xx_est = xx_est + (exp_prod * gamma[k, i])

def estimate_xx_norm(x_i, mu_k, sigma_k, m, o):
    assert x_i.ndim == 2 and mu_k.ndim == 2 and sigma_k.ndim == 2
    est_xx = sigma_k[np.ix_(m, m)] - np.dot(np.dot(sigma_k[np.ix_(m, o)], np.linalg.inv(sigma_k[np.ix_(o, o)])), np.transpose(sigma_k[np.ix_(m, o)]))
    est_xx = est_xx + np.dot(estimate_x_norm(x_i, mu_k, sigma_k, m, o), np.transpose(estimate_x_norm(x_i, mu_k, sigma_k, m, o)))
    assert est_xx.ndim == 2
    return est_xx

This work is based on this paper:

enter image description here

Related Questions

Updated April 18, 2018 10:19 AM

Updated June 06, 2017 00:19 AM

Updated July 09, 2018 11:19 AM

Updated February 15, 2019 12:19 PM

Updated July 03, 2017 15:19 PM