Support multidimensionality

[pavelkomarov]

In numpy we can often feed a method a multidimensional array and an axis parameter that specifies the direction we're operating along. The methods in this package are designed to work on univariate data streams, but it would be nice to generalize them to work with arrays of more dimensions.

Alternatively we could keep all methods univariate and force users to do something like apply_along_axis, however this loops up in Python, not taking advantage of speedups from numpy that may be possible with calls like fft, which push that loop down into compiled C.

[pavelkomarov]

This one will require going method by method, adding an axis parameter where possible. It's not terribly hard, but it will touch a lot of code. I recommend picking one and extending it to take an axis parameter and be capable of handling a multidimensional x. Then extend the test code to ensure it works. The first one will involve defining some new $f(x,y)$ in test_diff_methods to test against, but the others can reuse that same example.

[florisvb]

A challenge i see is that for many methods it may make sense to have different parameters for each dimension. This seems like something that is easy, and maybe better, for the user to handle. I would consider this low priority for now compared to many of the other issues.. Floris van Breugel | http://www.florisvanbreugel.com Assistant Professor of Mechanical Engineering & Graduate Program for Neuroscience University of Nevada, Reno Wildlife and Landscape Photography Galleries: http://www.ArtInNaturePhotography.com/ Blog: http://www.ArtInNaturePhotography.com/wordpress/

…

On Tue, Jul 29, 2025 at 3:30 PM Pavel Komarov ***@***.***> wrote: *pavelkomarov* left a comment (florisvb/PyNumDiff#76) <#76 (comment)> This one will require going method by method, adding an axis parameter where possible. It's not terribly hard, but it will touch a lot of code. I recommend picking one and extending it to take an axis parameters and optionally a multidimensional x, then extending test the test code to ensure it works. The first one will involve defining some new $f(x,y)$ in test_diff_methods to test against, but the others can reuse that same example. — Reply to this email directly, view it on GitHub <#76 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAB4EPHJO55LH5ILCEOCZ7T3K7YXPAVCNFSM6AAAAAB6D3J4DOVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZTCMZUGI2TEMBRGU> . You are receiving this because you are subscribed to this thread.Message ID: ***@***.***>

[pavelkomarov]

I envision you would still have to call one dimension at a time, specified by a scalar passed to axis, so you could separate that potential jumble by providing different parameters on each call.

[pavelkomarov]

@Jacob-Stevens-Haas, would you be interested in doing this one? I want to bring people in on this project and don't want to be the only one making changes, especially if the aim is to be compatible with another project you're more familiar with than I am.

[Jacob-Stevens-Haas]

I might have time to do this... How amenable would you be to a solution as in derivative? That is, temporarily reshape the data, apply differentiation along axis 0, iterating through an ndindex of the other axes, then revert the reshape? Pseudocode:

old_shape = temp_reshape(arr):
output = np.empty(arr.shape)
for multi_ind in np.ndindex(arr.shape[1:]):
    output[:, *multi_ind] = differentiate(arr[:, *multi_ind])
output = output.reshape(old_shape)

The considerations for this method are:

• Where would you want this happen here, e.g. a decorator to functions? in derivative all approaches are subclasses implemented by hooks, so it's easy to stick this kind of preprocessing into the user method in the superclass.
• Do you always return a same-shape array (no trimming endpoints)?
• I don't think runtime/memory efficiency matters much, since I assume most methods are dominated by large matrix operations or iterative optimization, but just something to be aware of

Testing is pretty easy: generate a random matrix, copy along one axis so that the derivative along that axis is always zero, verify it.

[pavelkomarov]

I find reshape calls confusing, personally. I'd rather index with [axis] if possible. It's a bit more direct. But if that can't work or is too complicated to figure out, reshaping is reasonable, because computationally it's essentially free.

Likewise, doing it with a decorator seems over-fancy to me if it can be done directly. A decorator would wrap all functions with the same layered-on loop, which to be common to everything probably would end up looking a lot like (and being spiritually redundant with) an apply_along_axis call. That wouldn't allow us to take advantage of any possible speedups that might come from pushing multidimensional data handling down into numpy. I agree those speedups aren't likely to be big in most cases, but I don't want to overgeneralize and rob us of the opportunity without turning some ingenuity on it and pondering first.

Yes, return values should always have the same dimension as the inputs.

I'd test with some intentional functions, like I do halfway down this notebook: https://github.com/pavelkomarov/spectral-derivatives/blob/main/notebooks/chebyshev.ipynb. We want to compare against easy-to-find analytic answers.

[Jacob-Stevens-Haas]

I forgot, yeah, apply_along_axis is right and probably what derivative should've used in the first place.

As for the decorator, two alternatives are manually adding apply_along_axis to each function, or add something like derivative's dxdt(method: str, axis: int, *args, **kwargs) that reads the method string and dispatches to the appropriate helper function, and can be reused for any other user-facing common guard code.

[pavelkomarov]

My instinct would be to give each method an axis: int=0 parameter individually rather than wrap with a common dxdt function. Unless it becomes clear that internally we are just calling along axis for all and there truly are no opportunities to do better, in which case a decorator or something that invokes apply_along_axis may make sense. I suppose you could decorate only the subset of methods where you can't do better.

I'd always prefer to pass a method directly by function handle than by name with a str.

How do docstrings work with decorators? Numpy functions typically have an axis parameter documented directly (e.g. https://numpy.org/doc/2.3/reference/generated/numpy.min.html), so to match that the docs should inherit and slot in a description. Going method by method gives finer-grained control over that, though I admit it's a lot of function headers to touch.

[Jacob-Stevens-Haas]

I'd always prefer to pass a method directly by function handle than by name with a str.

Users don't. solve_ivp, cvxpy.Problem, scipy.optimize.minimize all take a string for an integration method or optimization method, even though they also provide lower-level callables for each method. If I'm writing a BaseDifferentiation subclass, I don't want to make users import the correct callable, so I would have to all the possible imports myself, create a name for them, and add guard code. If a new method is added to pynumdiff it doesn't become available to pysindy until I patch imports.

OTOH, if pynumdiff does it, that process is seamless for users. derivative handles this with a decorator, @register(name=short_name) that adds a method to a global dictionary. Even if you just provide the dictionary (or an Enum), users could call derivative(arr, method_dict[shortname]).

How do docstrings work with decorators

Use functools.wraps