8
$\begingroup$

Non negative least square solves

$$ y = \alpha^T x \\ s.t. \forall i,~ \alpha_i \geq 0 $$

However, I would like to apply the non negativity constraint only on some coefficient, say only the $\alpha_i,~ \forall i \in \mathcal{S}$.

This is not possible using scipy's nnls function. However, I think I can solve the problem by using nnls with the following trick:

$$ y = \sum_{i \in \mathcal{S}} \alpha_i x_i + \sum_{i \notin \mathcal{S}} \beta_i x_i - \sum_{i \notin \mathcal{S}} \gamma_i x_i \\ s.t.~\forall i~, \alpha_i \geq 0, \beta_i \geq 0, \gamma_i \geq 0 $$

but this is hacky. Do you know any over way to implement this in python?

$\endgroup$

3 Answers 3

3
$\begingroup$

Your approach with the extra parameters will be problematic as multiple values for the coefficients will provide a same solution (the solution won't be unique).

Non-negative least squares regression is often solved by an active set method that updates in steps the active constraints (see for instance How do Lawson and Hanson solve the unconstrained least squares problem?). In those steps various regular least squares estimates are computed for different active sets. Those steps will fail with your approach because of the multiple solutions (you will get some errors during the computations like non-invertible matrices).

An alternative nnls for Python (aside from Galen's response) could be scipy.optimize.lsq_linear, which allows setting individual constraints and uses an active set method.

$\endgroup$
1
  • $\begingroup$ In the end, I prefered your solution as lsq_linear method of scipy was a clearer and most easy to use alternative. $\endgroup$
    – LucG
    Commented Nov 8, 2022 at 10:08
5
$\begingroup$

Here is an example of multiple linear regression in which some of the variables are unbounded while others are bounded. I'm going to have ten unbounded parameters and ten non-negative parameters in this example. The quick gist is to use scipy.optimize.Bounds in an optimizer that supports this argument such as 'trust-constr', along with the insight that it allows the use of np.inf. You just set some of the bounds to be -np.inf and np.inf for the unbounded parameters, and set the bounds to be 0 and np.inf for the non-negative parameters.

import numpy as np
from scipy.optimize import minimize, Bounds

np.random.seed(2022)

bounds = Bounds([-np.inf] * 10 + [0] * 10, [np.inf] * 20)
x = np.random.normal(size=100 * 20).reshape(100, 20)
true_a = np.arange(1,21)
true_y = x @ true_a + np.random.normal(size=100)
a0 = np.random.normal(size=20)

def f(a):
    y_hat = x @ a
    resid = true_y - y_hat
    lsq = np.mean(np.power(resid, 2))
    return lsq
    

result = minimize(f, x0=a0, method='trust-constr', bounds=bounds)

print(result.x)

Note that the last line accesses results.x rather than results.abecause x is priviliged in scipy.optimize's programming to be the parameter. But this prints a result for the optimized vector of parameters a as hoped. Here is the printout for this seed:

[ 0.9363189   2.01708136  3.07371156  4.03469035  4.97227273  5.91210627
  6.86926581  8.05433955  8.83633234  9.93401828 10.96973645 12.05863185
 12.95428506 13.9473809  14.93419422 16.05477142 17.13887755 18.37746539
 19.05598047 19.78565871]
$\endgroup$
13
  • 2
    $\begingroup$ It is one of several trust region methods available in scipy.optimize. The docs say their implementation is based on Lalee et al 1998 and Byrd et al 1997. I don't know which equivalence relation you are referring to, but I do not know any related results offhand anyway. $\endgroup$
    – Galen
    Commented Nov 4, 2022 at 18:00
  • 2
    $\begingroup$ At a glance in the papers I see gradients, so I suspect it is not derivative free, but you can read the papers I linked to check more carefully. $\endgroup$
    – Galen
    Commented Nov 4, 2022 at 18:00
  • 1
    $\begingroup$ This is one of the reasons that I dislike Scipy/Python, it is often difficult to find out what the underlying method is and the explanations of approaches are very superficial (little detail) with only references to articles that require a lot of research effort. In the end, this approach of using minimize is not so great because it approaches specific linear least squares problem as a general minimization problem. $\endgroup$ Commented Nov 4, 2022 at 18:48
  • 1
    $\begingroup$ SciPy has a GitHub where you can bring attention to issues and propose changes to the documentation. $\endgroup$
    – Galen
    Commented Nov 4, 2022 at 18:58
  • 1
    $\begingroup$ This gist shows that scipy.optimize.lsq_linear is very close in results to minimize(method="trust-constr"). $\endgroup$
    – Galen
    Commented Nov 5, 2022 at 4:22
1
$\begingroup$

The brute force approach would be using constrained minimization, applying the non-negativity constraints only to certain variables - along the lines of the python implementation provided in an answer to How to include constraint to Scipy NNLS function solution so that it sums to 1.

Another trick is variable transformation that makes certain variables always non-negative, like $x=y^2$. One can then use any general purpose minimization routine.

The only problem with these approaches is that they will perform much slower than NNLS for a large number of variables - whether this is a constraint depends on the problem, the available computational capacities, and the programming language used.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.