`first_value` investigation

[pavelkomarov]

The scipy code that performs integration returns an integration result one-shorter than the input, which is inconvenient. So in the utility code an extra value is prepended to that result vector to make the result same-length. It's a bit of a puzzle what should be added here, and the choice made by the existing code doesn't make perfect sense to me, because it chooses first_value = x[0] - dxdt_hat[0], and the units of x and its derivative don't match. It's possible this actually is a good choice, but while reviewing #77, @florisvb couldn't recall why this value was chosen, so I'd like us to do some experimentation or thinking and then justify this with a code comment, or change it.

[pavelkomarov]

I've been playing with this. dxdt = np.linspace(5,-75,6) results in a first_value of -8. If I plot the result, it looks like this:

Over it I've plotted the integral with the first value of 0, which looks like a more natural continuation to me, given the linearity of the derivative. But I suppose maybe the rationale here is that since $\frac{dx}{dt}$ is positive at its first value, we have to force in a little positive slope at that point. I'm just not sure this is the right way to account for this.

[pavelkomarov]

Okay, I think I have a principled idea. Here dxdt = [ 5., -11., -27., -43., -59., -75.]. If we let x = scipy.integrate.cumulative_trapezoid(dxdt), then this is [ -3., -22., -57., -108., -175.]. Let's FD that by x[1:] - x[:-1] = [-19., -35., -51., -67.].

These are intermediate values. -19 is halfway between -11 and -27, and -35 is halfway between -27 and -43. So by extrapolation we really want to have a value at the beginning of this thing that is halfway between 5 and -11.

The first_value that gives us back a value halfway between dxdt[0] and dxdt[1] when we [1:] - [:-1] is: x[0] - first_value = (dxdt[0] + dxdt[1])/2. That looks a little more like the code you had there before #77 (first_value = x[0] - np.mean(dxdt_hat[0:1])), although that code had an indexing-related bug (Matlab making an appearance) and is equivalent to the code that's currently there first_value = x[0] - dxdt_hat[0].

[pavelkomarov]

You know what's funny? That equation always results in first_value = 0, no matter the signal, because x[0] will already equal the average of dxdt[0] and dxdt[1], since that's how cumulative_trapezoid works with a default dx=1.0.

So intuition and math agree: The first value should be 0, and there's been a bug all this time!

[florisvb]

The units only dont match after dividing by dt, but I think at this step that hasnt been done yet? 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 Thu, Jun 5, 2025 at 4:12 PM Pavel Komarov ***@***.***> wrote: *pavelkomarov* left a comment (florisvb/PyNumDiff#88) <#88 (comment)> You know what's funny? That equation always results in first_value = 0, no matter the signal, because x[0] will already equal the average of dxdt[0] and dxdt[1], since that's how cumulative_trapezoid works with a default dx=1.0. So intuition and math agree: The first value should be 0, and there's been a bug all this time! — Reply to this email directly, view it on GitHub <#88 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAB4EPA5N74WR3IU3PHUIUD3CD2J5AVCNFSM6AAAAAB6UAIW6CVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDSNBXHAYTMNZVGY> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[pavelkomarov]

Some tests are no longer passing on #89, because I think they are just trying to match against values formerly generated by the code itself. I still believe a first value of 0 is actually correct, and that we need to correct all those tests. I've added test code for the integration utility grounded in my math above.

[pavelkomarov]

In _linear_model in the function __integrate_dxdt_hat_matrix__ there is another one of these motifs: first_value = x[:, 0:1] - np.mean(dxdt_hat[:, 0:1], axis=1).reshape(dxdt_hat.shape[0], 1).

Indexing [0:1] and taking the mean would give you the mean of the first two things in Matlab, but in Python this is a mistake. And analogous to the scalar case this issue is originally about, I believe the first_value should also be 0 here.

[florisvb]

I’m not positive, but its possible that the 0:1 syntax keeps the shape of the array, and maybe there are cases where that matters. As for being zero, I dont think that is correct. I think what this function is trying to do is extrapolate the first value based on the original first value and the derivative. I recall needing to check this very carefully by comparing the numerical derivatives with analytic derivatives, so if you make any changes please verify with analytic derivatives evaluated at the appropriate time stamps. 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 Fri, Jun 6, 2025 at 12:27 PM Pavel Komarov ***@***.***> wrote: *pavelkomarov* left a comment (florisvb/PyNumDiff#88) <#88 (comment)> In _linear_model in the function __integrate_dxdt_hat_matrix__ there is another one of these motifs: first_value = x[:, 0:1] - np.mean(dxdt_hat[:, 0:1], axis=1).reshape(dxdt_hat.shape[0], 1). Indexing [0:1] and taking the mean would give you the mean of the first two things in Matlab, but in Python this is a mistake. And analogous to the scalar case this issue is originally about, I believe the first_value should also be 0 here. — Reply to this email directly, view it on GitHub <#88 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AAB4EPAMW2NQ6RB3L2JYOBD3CIIT3AVCNFSM6AAAAAB6UAIW6CVHI2DSMVQWIX3LMV43OSLTON2WKQ3PNVWWK3TUHMZDSNJQHE3TOMZRGM> . You are receiving this because you were mentioned.Message ID: ***@***.***>

[pavelkomarov]

>>>dxdt_hat = np.ones((10,10))*np.arange(10)
array([[0., 1., 2., 3., 4., 5., 6., 7., 8., 9.],
       [0., 1., 2., 3., 4., 5., 6., 7., 8., 9.],
       [0., 1., 2., 3., 4., 5., 6., 7., 8., 9.],
       [0., 1., 2., 3., 4., 5., 6., 7., 8., 9.],
       [0., 1., 2., 3., 4., 5., 6., 7., 8., 9.],
       [0., 1., 2., 3., 4., 5., 6., 7., 8., 9.],
       [0., 1., 2., 3., 4., 5., 6., 7., 8., 9.],
       [0., 1., 2., 3., 4., 5., 6., 7., 8., 9.],
       [0., 1., 2., 3., 4., 5., 6., 7., 8., 9.],
       [0., 1., 2., 3., 4., 5., 6., 7., 8., 9.]])
>>>dxdt_hat[:,0:1]
array([[0.],
       [0.],
       [0.],
       [0.],
       [0.],
       [0.],
       [0.],
       [0.],
       [0.],
       [0.]])

You're right. That slicing does preserve the dimension. Then you go and take the mean along axis=1, which is the horizontal direction, so you again are taking the mean of single numbers. It's just the case from earlier but repeated for every row of a 2D matrix.

[pavelkomarov]

If I make that line first_value = x[:, 0:1] - np.mean(dxdt_hat[:, 0:2], axis=1).reshape(dxdt_hat.shape[0], 1), which takes the mean of the first two columns rather than the mean of a single column, and print first_value, then when I run the tests I see all zeros.

[florisvb]

It could be that the 1 should have been a 2, I don't know. It could also be that I tried 0:2 but that in the end 0:1, in which case the mean doesnt do anything, worked better in most cases. Again I think the real test is comparing the numerical result with analytical results.

Where does this utility function actually get used?

[pavelkomarov]

$ grep -IR . -e "__integrate_dxdt_hat_matrix__"
./pynumdiff/linear_model/_linear_model.py:def __integrate_dxdt_hat_matrix__(dxdt_hat, dt):
./pynumdiff/linear_model/_linear_model.py:    integral_X = __integrate_dxdt_hat_matrix__(X, dt)

So looks like it's only called in one place. The scalar version is called in more places.

I'll do some checks against analytic functions. #90 and #94 are meant to address that directly. We'll get to see whether any of the methods dependent on these integration functions end up with worsening or improving error bounds when we change first_value.