ENH: `linalg.solve_toeplitz`: add option to return reflection coefficients

[harrymander]

trafficstars

What does this implement/fix?

Exposes the reflection coefficients as returned by the internal levinson function. As explained in the docstring for levinson, this is useful for getting the partial autocorrelation function when estimating the parameters of an autoregressive process.

linalg.solve_toeplitz will return a 2-tuple of the solution and the reflection coefficients if the optional return_reflection=True is passed. The default is return_reflection=False, in which case the function behaves the same as at present.

[harrymander]

@j-bowhay Should I make a different function then (e.g. solve_toeplitz_with_reflection), so that the original function remains backwards-compatible?

[j-bowhay]

@j-bowhay Should I make a different function then (e.g. solve_toeplitz_with_reflection), so that the original function remains backwards-compatible?

Sorry, my original comment was made in error as this function originally only had one output so the approach I proposed will not work. I think your current approach will probably be fine however I will leave that to the judgement of one of the linalg regulars.

[mdhaber]

@harrymander were you still interested in this?

If so, @harrymander @j-bowhay do either of you know of another very simple method for solving these problems that we might want to have? If so, we could add a method argument to help us add additional output information. The difference is that:

• return_reflection is not an option that we really want to have, we don't want to permanently have the number of outputs depend on a parameter, and we don't want two cycles of deprecations to achieve the end goal;
• method may be an option that we actually want to have, and there is a single-cycle deprecation strategy you could use to get rid of the old behavior.

If the method is specified, solve_toeplitz would return a result object that can include additional information; if not (method=None), the old behavior would be preserved. However, in that case, you could immediately emit a warning against unspecified method. At the end of the warning period, you have some options, but in any case, you would no longer need to preserve the old behavior.

[harrymander]

@mdhaber I agree that changing the return type based on an argument is not ideal. However there seems to be some precedent for it already in scipy, e.g. optimize.curve_fit has a full_output argument that is by default False: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html, plus a few others: https://docs.scipy.org/doc/scipy-1.12.0/search.html?q=full_output. So full_output could be a better name for the argument.

As to whether the argument should be a bool or default None with a warning if not explicitly provided: I suppose that is up to you and other maintainers to decide if you would want to eventually deprecate the old return behaviour for a new major version release. Personally, I think the current behaviour is fine as a default, as the reflection coefficients are only needed for specific applications (e.g. AR processes).

IMO a better alternative would be to keep solve_toeplitz as is and expose a public levinson function that returns both the solution and the reflection coefficients since technically the reflection coefficients are an 'implementation detail' specific to the Levinson-Durbin recursion and not necessarily specific to solving a Toeplitz system in general (although the docs for solve_toeplitz do specify that the Levinson-Durbin algorithm is used).

[mdhaber]

However there seems to be some precedent

Yes! The problems we've had in these cases have confirmed the intuition that it should be avoided.

IMO a better alternative would be to... expose a public levinson function

Sure, that is the easiest route when there are backward compatibility concerns, but we also try to avoid having almost redundant functions. It's been about a year with no other comments, so there might not be enough support to warrant that. One possibility to avoid the redundancy would be to incorporate assume_a='toeplitz' as an option of solve (which we should probably do anyway). Then we could create levinson (or better name, if there is one) as you suggest and deprecate solve_toeplitz.