Default axis of differentiation is axis 1, rather than axis 0 or -1

[pavelkomarov]

In differentiation.py, the axis parameter defaults to 1 and is used to index things like X.shape[1]. I think it's more intuitive for the default direction of differentiation to be the first or last dimension, not the 2nd, especially for 3D or higher-dimensional data.

[Jacob-Stevens-Haas]

Also in the object method d and x

FWIW in pysindy we explicitly set the axis when calling derivative, so we don't have to worry about whether changing the default here will affect pysindy.

[andgoldschmidt]

Originally, derivative was written with "snapshot matrices" in mind, so X = [x_1 x_2 ... x_T] was the shape to use, and anything else was out of scope.

Perhaps this is too restrictive now. As long as downstream is safe, there is no reason to prefer the old way if a new approach enables additional features.

[pavelkomarov]

If I change the default to 0 in differentiation.py, the tests throw a bunch of errors from _align_axes. I'm not sure this really ever worked the way it should. Maybe some of the tests were implicitly relying upon the default. I'll investigate.

[pavelkomarov]

It seems like the assumption was we'd never work with greater than 2D data. Is this right?

[andgoldschmidt]

Yes, that's exactly right (see previous post about snapshot matrix assumption).

[Jacob-Stevens-Haas]

#41 added _align_axes to allow differentiating data that was greater than 2D. Looks like I omitted to add a docstring to that function, but it orients the desired axis as axis 1 and reshapes everything else into axis 0. There's some note about taking care with order of operations, so that the inverse operations in _restore_axes work correctly.

It's odd that changing the default in smooth_x and dxdt would cause a test failure, since this gets passed through Derivative.d and then the axis convention provided by _align_axes is what each method's implementation of compute() sees.

[pavelkomarov]

So even if the data is 3D or greater, we reshape to 2D, differentiate along axis 1, and shape back. I'm trying to picture how that works mathematically. I guess as long as the axis of interest is kept in tact, then all the vectors along it can be stacked into a big matrix, and then you can reshape back. This may or may not be easier than leaving it as is and trying to do whatever smoothing or differentiation operations along whichever dimension. The tradeoff is this extra reshaping code so that various methods can consistently expect axis 1.

I don't hate it. But I'm adding a docstring, and my preference would be to reshape such that differentiation always happens along axis 0.

[Jacob-Stevens-Haas]

how that works mathematically

Alternative view: what we consider an n-dimensional numpy array is really a contiguous list of memory locations along with information what size strides to take when indexing each of the n-dimensions. Reshape, moveaxis, and flatten don't really have great mathematical analogues. Moving an axis is just a swap of two elements in the metadata tuple. Flattening several axes is just forgetting some metadata, which is why _restore_axes needs it back in orig_shape.

My preference would be to reshape such that differentiation always happens along axis 0.

Changing the default value in the public API is easy; changing the backend convention, less so.