Minimizing scipy.stats.multivariate_normal.logpdf with respect to covariance

Question

I have a python script where I compute the value of a normal log-likelihood function for a sample of bivariate data using scipy's multivariate_normal.log_pdf. I am assuming the values of the sample means and variances, leaving only the sample correlation between the variables as the unknown,

from scipy.stats import multivariate_normal
from scipy.optimize import minimize

VAR_X = 0.4
VAR_Y = 0.32
MEAN_X = 1
MEAN_Y = 1.2

def log_likelihood_function(x, data):
    log_likelihood = 0
    sigma = [ [VAR_X, x[0]], [x[0], VAR_Y] ]
    mu = [ MEAN_X, MEAN_Y ]
    for point in data:
        log_likelihood += multivariate_normal.logpdf(x=point, mean=mu, cov=sigma)
    return log_likelihood

if __name__ == "__main__":
    some_data = [ [1.1, 2.0], [1.2, 1.9], [0.8, 0.2], [0.7, 1.3] ]
    
    guess = [ 0 ] 

    # maximize log-likelihood by minimizing the negative 
    likelihood = lambda x: (-1)*log_likelihood_function(x, some_data)
    
    result = minimize(fun = likelihood, x0 = guess, options = {'disp': True}, method="SLSQP")

    print(result)

No matter what I set as my guess, this script reliably throws a ValueError,

ValueError: the input matrix must be positive semidefinite

Now, the problem, by my estimation, seems to be scipy.optimize.minimize is guessing values that create a covariance matrix that is not positive definite. So I need a way to make sure the minimization algorithm throws away values that are outside the domain of the problem. I thought to add a constraint to the minimize call,

## make the determinant always positive
def positive_definite_constraint(x):
    return VAR_X*VAR_Y - x*x

Which is basically the Slyvester Criteron for the covariance matrix and would ensure the matrix is positive definite (since we know the variance is always positiv, that condition doesn't need checked) But it seems like scipy.optimize.minimize evaluates the objective function before it determines if the constraints are satisfied (which seems like a design flaw; wouldn't it be faster to search for a solution in a restricted domain, instead of searching all possible solutions and then determining if the constraints are satisfied? I might be mistaken about the order of evaluation, though.)

I am not sure how to proceed. I realize I am stretching the purpose of scipy.optimize here a bit by parameterizing the covariance matrix and then minimizing with respect to that parameterization, and I know there are better ways to calculate the correlation for a normal sample, but I am interested in this problem because of its generalization to distributions that are not normal.

Any suggestions? Is there a better way to solve this problem?

user_4745529 · Answer

You are on the right track. Note that your definiteness constraint reduces to a simple bound on the optimization variable, i.e. -∞ <= x[0] <= VAR_X*VAR_Y. Variable bounds are better handled internally than the more general constraints, so I'd recommend something like this:

bounds = [(None, VAR_X*VAR_Y)]
res = minimize(fun = likelihood, x0 = guess, bounds=bounds, options = {'disp': True}, method="SLSQP")

This gives me:

     fun: 6.610504611834715
     jac: array([-0.0063166])
 message: 'Optimization terminated successfully'
    nfev: 9
     nit: 4
    njev: 4
  status: 0
 success: True
       x: array([0.12090069])

Minimizing scipy.stats.multivariate_normal.logpdf with respect to covariance

1 Answers