3

I am trying to find weights across a number of forecasts to give a result that is as close as possible (say, mean squared error) to a known target.

Here is a simplified example showing three different types of forecast across four data points:

target = [1.0, 1.02, 1.01, 1.04]  # all approx 1.0
forecasts = [
    [0.9, 0.91, 0.92, 0.91],  # all approx 0.9
    [1.1, 1.11, 1.13, 1.11],  # all approx 1.1
    [1.21, 1.23, 1.21, 1.23]  # all approx 1.2
]

where one forecast is always approximately 0.9, one is always approximately 1.1 and one is always approximately 1.2.

I'd like an automated way of finding weights of approximately [0.5, 0.5, 0.0] for the three forecasts because averaging the first two forecasts and ignoring the third is very close to the target. Ideally the weights would be constrained to be non-negative and sum to 1.

I think I need to use some form of linear programming or quadratic programming to do this. I have installed the Python quadprog library, but I'm not sure how to translate this problem into the form that solvers like this require. Can anyone point me in the right direction?

Matthew Strawbridge
  • 19,940
  • 10
  • 72
  • 93
  • are there any restrictions on the coefficients (e.g., all positive, sum up to 1, etc)? How many vectors are there to combine, comparing to dimensionality? – Marat Sep 10 '21 at 20:11
  • Yes, the coefficients should all be non-negative. Ideally they should all sum to 1, although the ability to turn this constraint on and off would be even better. The length of each vector will exceed the number of vectors. – Matthew Strawbridge Sep 10 '21 at 20:43

2 Answers2

5

If I understand you correctly, you want to model some optimization problem and solve it. If you are interested in the general case (without any constraints), your problem seems pretty close to the regular least square error problem (which you can solve with scikit-learn for example).

I recommend to use cvxpy library for modeling an optimization problem. It's a convenient way to model a convex optimization problem, and you can choose which solver you want to work in the background.

Expanding cvxpy least square example, by adding the constraints you mentioned:

# Import packages.
import cvxpy as cp
import numpy as np

# Generate data.
m = 20
n = 15
np.random.seed(1)
A = np.random.randn(m, n)
b = np.random.randn(m)

# Define and solve the CVXPY problem.
x = cp.Variable(n)
cost = cp.sum_squares(A @ x - b)
prob = cp.Problem(cp.Minimize(cost), [x>=0, cp.sum(x)==1])
prob.solve()

# Print result.
print("\nThe optimal value is", prob.value)
print("The optimal x is")
print(x.value)
print("The norm of the residual is ", cp.norm(A @ x - b, p=2).value)

In this example, A (the matrix) is a matrix of all your vector, x (the variable) is the weights, and b is the known target.

EDIT: example with your data:

forecasts = np.array([
    [0.9, 0.91, 0.92, 0.91],
    [1.1, 1.11, 1.13, 1.11],
    [1.21, 1.23, 1.21, 1.23]
])

target = np.array([1.0, 1.02, 1.01, 1.04])
x = cp.Variable(forecasts.shape[0])
cost = cp.sum_squares(forecasts.T @ x - target)
prob = cp.Problem(cp.Minimize(cost), [x >= 0, cp.sum(x) == 1])
prob.solve()
print("\nThe optimal value is", prob.value)
print("The optimal x is")
print(x.value)

Output:

The optimal value is 0.0005306233766233817
The optimal x is
[ 6.52207792e-01 -1.45736370e-24  3.47792208e-01]

results are approximately [0.65, 0, 0.34] which is different from the [0.5, 0.5, 0.0] you mentioned, but that depends on how you define your problem. This is a solution for the least squares error.

Roim
  • 2,986
  • 2
  • 10
  • 25
  • Is there a reason you didn't use the data from the post i.e. target & forecasts to solve the example problem? – DarrylG Sep 10 '21 at 20:32
  • @DarrylG not really I just stayed with the example in the cvxpy documentation. I edited my answer and added example with the data from the post – Roim Sep 10 '21 at 20:53
2

We can see this problem as a least squares, which is indeed equivalent to quadratic programming. If I understand correctly, the weight vector you are looking for is a convex combination, so in least squares form the problem is:

minimize  || [w0 w1 w2] * forecasts - target ||^2
    s.t.  w0 >= 0, w1 >= 0, w2 >= 0
          w0 + w1 + w2 == 1

There is a least-squares function you can use out of the box in the qpsolvers package:

import numpy as np
from qpsolvers import solve_ls

target = np.array(target)
forecasts = np.array(forecasts)
w = solve_ls(forecasts.T, target, G=-np.eye(3), h=np.zeros(3), A=np.array([1, 1., 1]), b=np.array([1.]))

You can check in the documentation that the matrices G, h, A and b correspond to the problem above. Using quadprog as the backend solver, I get the following solution on my machine:

In [6]: w
Out[6]: array([6.52207792e-01, 9.94041282e-15, 3.47792208e-01])

In [7]: np.dot(w, forecasts)
Out[7]: array([1.00781558, 1.02129351, 1.02085974, 1.02129351])

Which is the same solution as in Roim's answer. (CVXPY is indeed a great way to start!)

Tastalian
  • 369
  • 4
  • 12