ENH: Add non-negative reconciliation heuristic for MinTraceSparse

[christophertitchen]

This PR adds a non-negative reconciliation heuristic for sparse minimum trace reconciliation, which currently is not supported in the fit method for MinTraceSparse. I have implemented something similar in my job as the time series we deal with are inherently non-negative in nature, as they are across many industries, and thought it would be useful. It is fast and guarantees non-negative coherence at the expense of producing an approximate solution with potentially ill-defined mathematical properties. However, the goal of such a transformation is to align forecasts with their use in practice for the users of our forecasts, not necessarily improve accuracy, which is a nice additional benefit that in my experience actually happens rather frequently.

We can expand this further in the future by implementing a sparse non-negative least squares solver—there seem to be a few nice examples of this in the signal processing space.

Step 1: Clip the base forecasts to remove any negative point forecasts, which of course does not guarantee non-negativity in the minimum trace method, but nonetheless provides a good starting point.

$\alpha_{h} = \begin{cases} \hat{y}_{i,h} & \text{if } \hat{y}_{i,h} \gt 0 \\ 0 & \text{otherwise} \end{cases}$

Step 2: Although it is now sufficient to ensure that all of the entries in $P$ are positive, as it is implemented as a linear operator, we solve the sparse linear system iteratively as usual using BiCGSTAB. We could also try changing this to MINRES to exploit the symmetry of $\left( S'W^{-1}_{h}S \right)^{-1}$ in the future which has been more performant in my experience, but it can be finicky to find the relative tolerance that best balances speed with maximum error.

$\beta_{h} = S\left( S'W^{-1}_{h}S \right)^{-1}S'W^{-1}_{h}\alpha_{h}$

Step 3: We now have reconciled point forecasts which minimise the error variances of the coherent forecasts, but which may be negative, so we need to check if any of the coherent bottom level point forecasts are negative. If so, the bottom level negative forecasts are clipped.

$b_{h} = \begin{cases} \beta_{j,h} & \text{if } \beta_{j,h} \gt 0 \\ 0 & \text{otherwise} \end{cases}$

$j = \left\{ n_{a}, n_{a} + 1, \cdots , n_{a} + n_{b} - 1, n_{a} + n_{b} \right\}$

Step 4: Restore coherence by aggregating the non-negative bottom level reconciled forecasts across the data structure, which leads to an approximate solution as we used the leverages of the negative-containing $SP$ in the minimum variance unbiased estimator to generate $b_{h}$.

$\tilde{y}_{h}= Sb_{h}$

Step 5: This directly forces projection onto the non-negative coherent subspace, which can be equivalently and subsequently achieved by using the bottom-up reconciliation projection matrix. We also set $W$ appropriately for probabilistic reconciliation if required.

$\tilde{y}_{h} = S\begin{bsmallmatrix} 0_{n_{b}, n_{a}} \quad I_{n_{b}} \end{bsmallmatrix} \tilde{y}_{h}$

$W = I_{n}$

[review-notebook-app]

Check out this pull request on 

See visual diffs & provide feedback on Jupyter Notebooks.

---

Powered by ReviewNB

[elephaint]

Great work again - I took the liberty to review it. Minor comment in the code.

Apart from that, I'm sure you checked it but it may be worth checking the types of the sparse matrices in the various matrix multiplications (some may be csc, others csr) and whether it's better to transpose operations or not. This can often lead to easy speedups.

Other than that, I think mostly what is needed are a few unit tests, comparable to the non-sparse MinT tests.

[elephaint]

Do you remove this as SciPy no longer supports it? (I don't see it being mentioned as a valid value on bicgstab page)

[elephaint]

There was an issue on Github with this too, I think this argument is indeed deprecated and needs to removed.