constrained convex optimisation for a multivariate normal distribution? Given some multivariate normal distribution $P(\mathbf{x|\mu,\Sigma} )$ I'm looking for a method to maximise $P$ subject to $x_i \ge 0 \; \forall \; i$.
I guess this is typically approached by instead minimising
$f(\mathbf{x}) = (\mathbf{x - \mu})^{\top}\mathbf{\Sigma}^{-1}(\mathbf{x - \mu}) $
subject to our constraints on $\mathbf{x}$.
Can anyone suggest an algorithm for solving the problem?
I'm happy to implement it myself, since I need it to be in Python and use the scipy.sparse matrix formats (as $\mathbf{\Sigma}$ is large and sparse).
Thanks!
 A: The optimization problem$$\min_x (x-\mu)^T\Sigma^{-1} (x-\mu)$$
subject to $x \geq 0$ falls under a class of problem known as quadratic programming problem.
A: Why not just use lognormals instead of normals if you want to force x>0? Then you can just solve it with maximum likelihood estimation.
Numerically, the python package CVXOPT should be able to do this and supports sparse matrices. You can figure out how to change matrix to sparse matrix in a way that suits your problem here:
http://cvxopt.org/userguide/matrices.html
Example code for dense matrices:
def opt(mu, cov):
    """
    minimize    (x -mu )' Q (x - mu)
    subject to  x > 0
    where Q is the precision matrix cov^(-1)

    CVXOPT minimizes    x'Px + q'x·
    subject to          Gx <= h
                        Ax == b

    The objective can be rewritten as x'Qx - 2 mu' Q x
    since mu' Q mu is a constant.

    So gathering up the terms we have that P = Q,·
    q' = -2 mu' Q, G = -1, A = 0, b = 0 and h = 0

    mu: means of shape n_samples
    cov: covariance matrix of shape (n_features, n_features)

    """

    n_features = cov.shape[0]
    Q = np.linalg.inv(cov)

    # objective
    P = matrix(Q)
    q = matrix((-2*mu[None, :].dot(Q)).T)

    # constraints

    G = matrix(-np.identity(n_features))
    A = matrix(0.0, (n_features, n_features))
    b = matrix(0.0, (n_features, 1))
    h = matrix(0.0, (n_features, 1))

    sol = solvers.qp(P, q, G, h, A, b)
    return np.asarray(sol['x']).ravel()

However this gives rank errors for some reason, meaning that something with the problem is ill defined (or that I made some error). Hope it helps at least.
