I am new to constraint optimization and I am wondering if it is possible to do a constraint maximum likelihood estimation for a least square problem with cvxpy?

My basic idea is minimize the sum of squares by fitting N distributions (the type of distribution does not matter but ideally similiar to Gaussian distributions).

So, I was thinking I could define the distribution parameters as cp.variables and use the cp.variables and a given domain_1xP to calculate for each distribution the corresponding pdf values W_1xP. Then I could stack them and use them as a matrix to minimize my problem.

However, this seems to be not possible due to Disciplined Convex Programming (DCP) requirements. I was wondering if there is some kind of a workarround to define a DCP compliant distributions based on cp.variable parameters.

I tried something like this:

# Define distribution parameters as cp variables
# N = number of technologies, G = number of Gaussians per technology
means_NxG = cp.Variable((N, G))
stddevs_NxG = cp.Variable((N, G), nonneg=True)
scaling_factors_NxG = cp.Variable((N, G), nonneg=True)

# Calculate W_NxP based on distribution parameters
# P = number of points in the price domain
W_NxP = calculate_with_G_normal_distributions_W_NxP_cvxpy(
    price_domain_Px1=price_domain_Px1, 
    means_NxG=means_NxG, 
    stddevs_NxG=stddevs_NxG, 
    scaling_factors_NxG=scaling_factors_NxG
)

# Minimize sum of squares
objective_TxP = y_TxP - ACF_TxN @ W_NxP
objective_loss_1x1 = cp.norm(objective_TxP, 'fro')
objective = cp.Minimize(objective_loss_1x1)

prob = cp.Problem(
    objective=objective, 
    constraints=[some constraints including W_NxP]
)
prob.solve(verbose=True)

And call the function:

def calculate_with_G_normal_distributions_W_NxP_cvxpy(
        price_domain_Px1: Px1,
        means_NxG: cp.Variable, 
        stddevs_NxG: cp.Variable, 
        scaling_factors_NxG: cp.Variable
):
    N, G = means_NxG.shape  
    P = price_domain_Px1.shape[0]  # P = number of points in the price domain
    
    # Initialize W_NxP as a list that will hold the results for each technology
    W_NxP = []

    # Flatten the price domain
    price_domain_Px_ = price_domain_Px1.flatten()

    # Build W_NxP using cvxpy expressions
    for technology in range(N):  # Loop over technologies
        technology_density_1xP = 0  # Start with zero for each technology
        for gaussian_parameter in range(G):  # Loop over Gaussian components for each technology
            mean = means_NxG[technology, gaussian_parameter]
            stddev = stddevs_NxG[technology, gaussian_parameter]
            scaling = scaling_factors_NxG[technology, gaussian_parameter]

            # Calculate the Gaussian pdf for each point in the price domain
            pdf_values_1xP = scaling * cp.exp(-0.5 * cp.square((price_domain_Px_ - mean) / stddev))

            # Sum up the pdf values for each Gaussian in this technology
            technology_density_1xP = technology_density_1xP + pdf_values_1xP

        # Reshape the result into a 1xP row vector and append it to the list
        W_NxP.append(cp.reshape(technology_density_1xP, (1, P)))

    # Stack the list of 1xP vectors vertically to create the final NxP matrix
    W_NxP = cp.vstack(W_NxP)
    return W_NxP

I know that there are different solutions and I implemented a working version using gradient descent but I like the non parametrized cvxpy solutions, e.g. if I define a W_NxP = cp.Variable(N,P) it works very well and the output looks a combination of normal distributions. SO, I was wondering if there is a way to make this DCP compliant?