Non negative least square on some coefficient

Question

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?

Answer 1

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.

Answer 2

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]

Answer 3

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.