I am using fmin_slsqp to find the weights that minimize mean squared error. The weights need to be positive. For each pair of X and y, it takes ~10 seconds. (Each X is (10, 1000) and y is (10,)). I have 8000 pairs that need to be calculated Is there any error with the code, or it is just my data takes too long to converge? Is there any way to make this process efficient, like is there a way to calculate all 8000 pairs together? def loss(W, X, y): return np.mean((y - X.dot(W))**2) def get_result(X, y): w_start = [1/X.shape[1]] * X.shape[1] weights = fmin_slsqp(partial(loss_w, X=X, y=y), np.array(w_start), bounds=[(0.0, np.inf)] * X.shape[1], disp=False) return weights Continue reading...