The Problem: Prediction Combination

Imagine you have several different models producing predictions on a regression task (temperature forecasting, traffic, bitcoin price, etc.). To get a final prediction, a simple way to combine the predictions would be a linear combination

where the weights of each forecast are determined by a simple ordinary least squares (OLS) fit using historical data. You could of course use more complicated modeling methods like decision trees or neural networks, but OLS is a simple and, in some circumstances, interpretable method. It also has the advantage of being easy to refit if one of the models’ predictions is missing.¹

Instability of OLS

However, fitting a simple OLS can run into problems. Consider the simple situation where the true underlying signal is

and there are two different models, the first of which predicts

and the second predicts

Notice that both predictions are highly correlated with the underlying target signal

If we fit an OLS using these two models, we get

Even though both predictions are positively correlated with the true target, the OLS gives negative weight to the first prediction. Intuitively, this feels wrong. Why would we take a useful prediction and then use its opposite direction to predict the target?

Geometry of Instability

One way to see what’s going on is to generate \(n\) points sampled according to the true distribution. Then compute the two different models’ predictions in order to form the \((n, 2)\) OLS problem. We can repeat this experiment many times, say 500, and plot the different weights resulting from the OLS as a scatter plot, together with the non-negative orthant.² Figure 1 below does exactly this.

Notice that for small \(n\), the variance in the weights is very high. This means that slight perturbations in the data lead to wildly different coefficients. This sensitivity has real consequences. Data is never measured to arbitrary precision: temperature rounded to two decimal places introduces rounding errors that can shift near-collinear OLS weights into a completely different region.

Rolling refits face the same problem. Swapping one day into a 500-day lookback window is a 0.4% change to the dataset, yet the widget shows the weights can change drastically and even change sign between consecutive refits.

The fundamental reason for the instability is that the two models’ predictions are highly correlated:

Stacking each model’s \(n\) predictions as the columns of a matrix \(A\), the minimizer of \(\frac{1}{2}\|A\beta-y\|_2^2\) is given by \(\beta^\star = (A^TA)^{-1}A^Ty\). \(A^TA\) is the Gram matrix of the predictor columns. In this particular case, the population Gram matrix is

with eigenvalues \(\lambda_1 \approx 14.93\) and \(\lambda_2 \approx 0.067\), a condition number of \(\kappa \approx 223\). This essentially says that the columns \(f_1\) and \(f_2\) are nearly collinear.

Figure 2 below makes this collinearity visible, plotting \(f_1\) and \(f_2\) as nearly-parallel arrows in feature space.

NNLS

Instead of solving the ordinary least squares problem, we can add a constraint that we expect the combination coefficients to be non-negative. This is Non-Negative Least Squares (NNLS).

As before, \(A \in \mathbf{R}^{n \times K}\) has one column per forecast, so \(\beta \in \mathbf{R}^K\) is the vector of combination weights and \(y \in \mathbf{R}^n\) is the target.

Put another way, we want to find the conic combination

that is closest to \(y\) (in the least squares sense) where each weight \(\beta_i\ge 0\).³

This feasible set is exactly the shaded wedge in Figure 2. The wedge is the cone of all non-negative combinations \(\beta_1 f_1 + \beta_2 f_2\) reachable under this constraint.⁴ Since \(y\) lies just outside that cone, no choice of non-negative weights can reconstruct it exactly. The best we can do is find the point in the cone closest to \(y\), which is precisely what NNLS solves for.

Figure 3 plots the level curves of the least squares objective function. Since the minimum sits outside the non-negative orthant, NNLS instead re-minimizes the objective along the boundary \(\beta_1 = 0,\ \beta_2 \geq 0\), landing at \((0, 8/13) \approx (0, 0.62)\).

It’s clear that our new NNLS solution achieves a worse objective value than OLS. Even though having non-negative weights when combining predictions that are correlated with the target makes intuitive sense, what benefit does it actually confer?

Figure 4 replays the same 500 repeated samples as Figure 1, but solves the NNLS problem each time.

We can see from the figure that the variance of the weights is massively reduced. In practice, this means small perturbations in data (e.g. numerical precision, rolling a sliding window of history forward) will not drastically change the computed weights. But how do we find the NNLS solution?

Solving NNLS

A natural first instinct is to solve the unconstrained OLS problem and clamp any negative weights to zero.

Figure 5 shows the OLS solution, the NNLS solution, and the clamped OLS solution.

From the figure we can see that clamping the OLS solution from \((-1,1)\) to \((0,1)\) leads to a sub-optimal point compared to \((0,8/13)\).

We can also compare NNLS vs clamping by viewing the predictions as vectors in \((X_1, X_2)\) space (as opposed to \((f_1,f_2)\) space), shown in Figure 6. On the left, OLS reaches \(y\) exactly by walking \(-f_1\) then \(f_2\) tip to tail, a path that leaves the cone entirely. On the right, clamping \(\beta_1\) to zero simply deletes the \(-f_1\) leg, landing at the tip of \(f_2\) instead of at \(y\). NNLS, on the other hand, finds the closest point to \(y\) that still lies in the cone.

Projected Gradient Descent

The simplest NNLS algorithm alternates two steps:

1. Gradient step: move in the steepest-descent direction

   \[\begin{align*} \tilde\beta &= \beta_t - \eta \nabla L(\beta_t) \\ &= \beta_t - \eta\,A^T(A\,\beta_t - y) \end{align*}\]

2. Project: clamp any negative weight back to zero

   \[\beta_{t+1} = \max(\tilde\beta,\ 0)\]

Figure 7 runs both the unconstrained and projected versions from the same starting point \(\beta_0 = (0, 0)\), iterating until a step changes \(\beta\) by less than \(10^{-5}\).

Lawson-Hanson Method

Another approach is to try to find directly which weights need to be positive and which need to be exactly zero. We can partition the weights into

• Passive set \(\mathcal{P}\): weights free to be positive
• Active set \(\mathcal{Z}\): weights clamped to zero (i.e. the constraint \(\beta_i \ge 0\) is active)

The first step of the algorithm is to initialize \(\mathcal{P} = \emptyset\) and \(\mathcal{Z} = \{1,\ldots,K\}\). Put another way, all \(\beta_i\) are initialized to 0.

In the outer loop of the algorithm we compute

which gives the vector of sensitivities of the loss with respect to each \(\beta_i\). Since we are trying to minimize the loss, it is more natural to work with the negative gradient

If \(z_i \le 0\) for every \(i \in \mathcal{Z}\), then we terminate since decreasing the loss would require making a \(\beta_i\) negative. Otherwise, we take the index of the largest component of \(z\) in the active set

and move it to the passive set \(\mathcal{P}\). Let \(A_{\mathcal{P}}\) be the submatrix of \(A\) containing only the columns indexed by \(\mathcal{P}\). We then solve the OLS problem

If the coefficients from the unconstrained OLS are all positive, we simply move to the next iteration and choose the next variable to add with the largest \(z_i\).⁵

Figure 8 shows the first few steps of the Lawson-Hanson algorithm on the 3-variable problem

The final step of Figure 8 shows that moving a variable into the passive set can cause the OLS solve to return a negative weight for a variable already in \(\mathcal{P}\).

We can perform a correction to make the parameters feasible again. Call the infeasible OLS proposal \(s\). For each \(i \in \mathcal{P}\) with \(s_i \le 0\), find the scalar \(\alpha_i\) such that

Geometrically, this is finding the scalar such that walking from the current feasible \(\beta\) toward the infeasible proposal \(s\) sets coordinate \(i\) to 0.

We can solve this equation analytically for each \(\alpha_i\):

Using \(\alpha_{\min} = \min_i \alpha_i\) ensures

This then becomes the new \(\beta\) (which is feasible by construction) and the coordinate corresponding to \(\arg \min_i \alpha_i\) becomes exactly zero so we move it from \(\mathcal{P}\) to \(\mathcal{Z}\).

Figure 9 lets you drag \(\beta\) (kept feasible) and the infeasible proposal \(s\) around the plane and observe the resulting \(\alpha_{\min}\).

With this correction for infeasible proposals in hand, Figure 10 continues on from Figure 8, showing the remainder of the algorithm.

We can be sure that we’ve reached the minimum when⁶

1. every \(\beta_i\) in the active set \(\mathcal{Z}\) has \(z_i \le 0\)
2. every \(\beta_i\) in the passive set \(\mathcal{P}\) has \(z_i = 0\)

If \(\beta_i\) were in \(\mathcal{Z}\) at optimality and had \(z_i > 0\), then we could decrease the loss by setting \(\beta_i=\delta\) (for sufficiently small \(\delta\)), which contradicts our original weight \(\beta_i\) being optimal.

Similarly, if \(\beta_i > 0\) were in \(\mathcal{P}\) at optimality and had \(z_i \ne 0\), we could move \(\beta_i\) by a tiny amount \(\delta\) in the direction of \(z_i\) without violating the constraint \(\beta_i \ge 0\). This again violates our assumption that \(\beta_i\) was optimal to begin with.

Figure 10’s final (done) step satisfies both of these conditions since every weight remaining in \(\mathcal{Z}\) has \(z_i \le 0\).

Conclusion

Using highly correlated features leads the OLS problem to be ill-conditioned. This manifests as high variance in the solution weights. Non-negative least squares effectively reduces the variance in exchange for a more biased solution.

While clamping the OLS weights to 0 is an intuitive way to get a weight vector with only non-negative weights, it turns out to be sub-optimal. Both projected gradient descent and the Lawson-Hanson method are algorithms for correctly solving the NNLS problem.