The best way to scale matrix variables in a SCIPY linear programming scheme

Question

I have the following optimization scheme implemented in NNLS in scipy .

import numpy as np from scipy.optimize import nnls from scipy import stats #Define problem A = np.array([[60., 90., 120.], [30., 120., 90.]]) b = np.array([6700.5, 699.,]) # Add ones to ensure the solution sums to 1 b = np.hstack([b,1.0]) A = np.vstack([A,np.ones(3)]) x, rnorm = nnls(A,b) print x # the solution is # [ 93.97933792 0. 0. ] # we expect it to sum to 1 if it not skewed 

As you can see, the vector b much higher than the values ​​in A My question is the best / reasonable way to scale A and b , so the solution is not skewed.

Note that both A and b are raw gene expression data without pre-processing.

Answer 1

If you want to enable the equality constraint, you cannot really use the nnls procedure, since it does not satisfy the equalities. If you are limited to what is offered in scipy, you can use this:

 import numpy as np from scipy.optimize import minimize #Define problem A = np.array([[60., 90., 120.], [30., 120., 90.]]) b = np.array([6700.5, 699.,]) #----------------------------- # I tried rescaling the data by adding this two lines, # so that they're in same scale. # but why the solution is different? # x: array([ 1., 0., 0.]) # What the correct way to go? #----------------------------- # A = A/np.linalg.norm(A,axis=0) # b = b/np.linalg.norm(b) def f(x): return np.linalg.norm(A.dot(x) - b) cons ={'type': 'eq', 'fun': lambda x: sum(x) - 1} x0 = [1, 0, 0] # initial guess minimize(f, x0, method='SLSQP', bounds=((0, np.inf),)*3, constraints=cons) 

Output:

  status: 0 success: True njev: 2 nfev: 10 fun: 6608.620222860367 x: array([ 0., 0., 1.]) message: 'Optimization terminated successfully.' jac: array([ -62.50927734, -100.675354 , -127.78314209, 0. ]) nit: 2 

This minimizes the objective function directly, and also creates a limit on the equality that interests you.

If speed is important, you can add information about jacobian and hessian or, even better, use the correct QP solver provided by cvxopt .