Hi, I am reaching out for help. First, the mathematical problem I am trying to solve: Suppose there is a rectangular matrix M = [[m11, m12, m13, m14], [m21, m22, m23, m24], [m31, m32, m33, m34]], from which I only know the totals (i.e., the sum of elements) by row and column. Hence, suppose my data is: tot_of_rows = [80.0, 230.0, 132.0] tot_of_cols = [74.0, 200.0, 91.0, 77.0] Problem: Infer the elements of the matrix using the method of Maximum Entropy. (See Cho, W. and Judge, G. (2008), "Recovering vote choice from incomplete data". Journal of Data Science 6, 155-171) Now, my current code that, in my mind, should solve this problem is (see also http://nbviewer.ipython.org/gist/anonymous/eee0f44d3de1a09570b0): # ----------------- beginning of python code ---------------------- # Importing libraries import numpy as np from scipy.optimize import minimize # Defining functions def Entropy(p_vec): # Function receives a vector of conditional probabilities, that come # from flattening a matrix of conditional probabilities. return -1.0 * np.dot(p_vec, np.log(p_vec)) def vec2mat(vec, nrows, ncols): # Re-generating the matrix from the vector mat = np.array([ [vec[j + i*ncols] for j in range(ncols)] for i in range(nrows) ]) return mat def mat2vec(mat): # Flattening the matrix into a vector nrows = int(np.array(mat).shape[0]) ncols = int(np.array(mat).shape[1]) vec = [mat[i, j] for i in range(nrows) for j in range(ncols)] return np.array(vec) # Functions for the Constraints: # Probabilities sum up to one # Probability: P[p,c] = "conditional probability of country c exporting p" # Xcp[c,p] = P[p,c]*Xc[c] # The problem is to maximize Entropy(P) # Subject to: # 1) sum_p P[p,c] - 1 = 0, for all c # 2) sum_c P[p,c]*Xc[c] - Xp[p] = 0, for all p # 3) 0 <= P[p,c] <= 1, for all c, p # for each c def Constraint_Norm(p_vec, Xc_vec, Xp_vec, c): lenvec_p = len(Xp_vec) lenvec_c = len(Xc_vec) # For clarity, reconstruct the matrix from p_vec p_mat = vec2mat(p_vec, lenvec_p, lenvec_c) return( np.sum([ p_mat[p,c] for p in range(lenvec_p) ]) - 1 ) # About the totals # for each p def Constraint_Mean(p_vec, Xc_vec, Xp_vec, p): lenvec_p = len(Xp_vec) lenvec_c = len(Xc_vec) # For clarity, reconstruct the matrix from p_vec p_mat = vec2mat(p_vec, lenvec_p, lenvec_c) return( np.sum([ p_mat[p,c]*Xc_vec[c] for c in range(lenvec_c) ]) - Xp_vec[p] ) # Initializing the parameters and data of the problem Xc = np.array([80.0, 230.0, 132.0]) Xp = np.array([74.0, 200.0, 91.0, 77.0]) size_c = len(Xc) size_p = len(Xp) # Initializing the tuple of constraints cons = () # Including constraint 1 cons += tuple( ({"type": 'eq', 'fun': lambda x: Constraint_Norm(x, Xc, Xp, c)} for c in range(size_c)) ) # Including constraint 2 cons += tuple( ({"type": 'eq', 'fun': lambda x: Constraint_Mean(x, Xc, Xp, p)} for p in range(size_p)) ) # Generating the bounds for each probability zeroonebounds = tuple((0,1) for c in range(size_c) for p in range(size_p)) # Generating the initial vector to start the 'minimize' algorithm init_vec = np.array([ 1.0/size_p + (np.random.random()-1)/10 for p in range(size_p) for c in range(size_c) ]) ##################################### # Finally, running the algorithm result = minimize(lambda x: -1.0 * Entropy(x), init_vec, bounds=zeroonebounds, constraints=cons, options={'disp': True}) # Display results res_mat = vec2mat(result.x, size_p, size_c) print res_mat # ----------------- end of python code ---------------------- This throws the following error: "Singular matrix C in LSQ subproblem (Exit mode 6)" Why? I don't see anything that would make a matrix in this context non invertible. Any help is highly appreciated.