Best options for SINDYc

[anfh27]

I get a perfect fit on my derivative using SindyC. However, when I simulate forwards, the prediction is not good. Changing regularization doesn't seem to have an impact here. Reducing the threshold does help to regularise but the predictions are still too flat. I do need to simulate forwards.

    diff =  ps.SmoothedFiniteDifference(drop_endpoints=True)
    opt = ps.SR3(threshold=0.00001, nu=100000000, max_iter=1500)
    fl = ps.PolynomialLibrary(degree = 1, include_bias=(True), include_interaction = True)
    model = ps.SINDy(optimizer=opt, feature_library=fl, differentiation_method=diff)

difference between model.differentiate and model.predict:

simulation:

It's easier to see what's happening if I include more components of the PCA:

Any idea what's going wrong?

[Jacob-Stevens-Haas]

Any idea what's going wrong?

Yeah, you didn't label your plots

[anfh27]

The control inputs are quite varied, their coefficients are not 0, and the derivative has a good fit.
I don't understand why the predictions are so bad.

[Jacob-Stevens-Haas]

To amplify my above answer:
first plot: Why would differentiate and predict be related? Also look at the scale: 1e-17. That means they're both approximately zero.
second plot: If the noisy lines are measured test data and the smooth lines are simulations starting at zero, this is a very good fit. Again - label your plot
third plot: I still don't know what is what, so label your plot

[anfh27]

first plot: I understood that model.differentiate calculates the derivative from the state, model.predict predicts the derivative based on pysindy output, so if the difference between them is zero, it's a good fit on the derivative.

second plot: Yes, the noisy lines are the measured data and the smooth lines are the simulations. The training phase ends at around 90 there, and the data fit with future control inputs only from ~90 onwards. The training data fits very well but the test phase is poor.

third plot: Not particularly important, it's just showing that the test phase predictions (again beginning from around 90) are straight lines instead of continuations of the same kind of dynamics as the training phase. I tested a few different systems and got similar results, but can't figure out the issue.

Apologies for the confusion.

[Jacob-Stevens-Haas]

First plot

Are the two series predict and differentiate? Or are they predict-differentiate for different coordinates?

second & 3rd plot

It doesn't look reasonable that the simulations were so close to the training data. Are you sure that's what you're plotting. Also, it looks like you don't have any test data? You are correct "straight lines instead of continuations of the same kind of dynamics as the training phase" is odd. I would suspect user error, but you didn't share any code for how the series are generated, and you didn't share model.print().

[anfh27]

I reproduced the error on a simple system. It seems like the problem happens when using an SVD matrix in place of controls. That's unfortunate, because I'm using a few delays on my real controls and it reliably produces a model instead of a "matrix is not positive definite" error.
I tried using StableLinearSR3, this didn't help.
Do you know any way of doing this?

Here's a test system:

import numpy as np
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import numpy.linalg as npl
import sklearn
import numpy.random as rnd
np.random.seed(42)
from sklearn.decomposition import PCA
import pysindy as ps


length = 130
num_variables = 4

random_variables = np.zeros((num_variables, length))

slopes = np.linspace(-0.01, 0.01, num_variables)


for i in range(num_variables):
    smooth_curve = np.sin(np.linspace(0, 2 * np.pi, length)) * np.random.uniform(0.5, 1.5)
    trend = slopes[i] * np.arange(length)
    random_variables[i, :] = smooth_curve + trend


for i in range(num_variables):
    complex_curve = (np.sin(np.linspace(0, (2 + i) * np.pi, length)) * np.random.uniform(0.5, 1.5) +
                     np.cos(np.linspace(0, (4 + i) * np.pi, length)) * np.random.uniform(0.3, 0.8) +
                     np.sin(np.linspace(0, (6 + i) * np.pi, length)) * np.random.uniform(0.2, 0.6))
    
    trend = slopes[i] * np.arange(length)
    random_variables[i, :] = complex_curve + trend

plt.figure(figsize=(10, 6))
for i in range(num_variables):
    plt.plot(random_variables[i, :], label=f'Variable {i+1}')
plt.title("Example System")
plt.xlabel("Index")
plt.ylabel("Value")
plt.legend()
plt.show()

X1 = random_variables.T
C1 = X1[1:]
pred = 8
timestep = 90
Xa = np.array(X1[:timestep,:])
u = C1[:timestep+pred]

diff = ps.SmoothedFiniteDifference()
opt = ps.SR3(threshold=0.001, nu = 100000, max_iter=1000)
fl = ps.PolynomialLibrary(degree=1, include_bias=True, include_interaction=True)
model = ps.SINDy(optimizer=opt, feature_library=fl, differentiation_method=diff)
t = np.linspace(0.0, len(Xa)-1, num=len(Xa))
t2 = np.linspace(0.0, len(t)+pred-1, num=len(t)+pred)           
model.fit(Xa, t, u=u[:len(Xa),:])

model.print()
u0 = Xa[0,:]
sim = model.simulate(u0, t2, u=u[:,:])

for i in range(num_variables):
    plt.plot(sim[:, i], label=f'Simulated Variable {i+1}')
    plt.plot(X1[:len(sim), i], '--', label=f'Original Variable {i+1}')

plt.title("SINDy Simulated vs. Original Data")
plt.xlabel("Time Index")
plt.ylabel("Variable Value")
plt.legend()
plt.show()

###### With U matrix of SVD as controls
u,s,v = npl.svd(C1[:timestep+pred], full_matrices=1)
opt = ps.SR3(threshold=0.001, nu = 100000, max_iter=1000)
fl = ps.PolynomialLibrary(degree=1, include_bias=True, include_interaction=True)
model = ps.SINDy(optimizer=opt, feature_library=fl, differentiation_method=diff)
t = np.linspace(0.0, len(Xa)-1, num=len(Xa))
t2 = np.linspace(0.0, len(t)+pred-1, num=len(t)+pred)           
model.fit(Xa, t, u=u[:len(Xa),:])

model.print()
u0 = Xa[0,:]
sim = model.simulate(u0, t2, u=u[:,:])

for i in range(num_variables):
    plt.plot(sim[:, i], label=f'Simulated Variable {i+1}')
    plt.plot(X1[:len(sim), i], '--', label=f'Original Variable {i+1}')

plt.title("SINDy Simulated vs. Original Data")
plt.xlabel("Time Index")
plt.ylabel("Variable Value")
plt.legend()
plt.show()


#### With Stable Linear SR3
opt = ps.StableLinearSR3(threshold=0.001, nu = 100000, max_iter=1000, verbose=True)
fl = ps.PolynomialLibrary(degree=1, include_bias=True, include_interaction=True)
model = ps.SINDy(optimizer=opt, feature_library=fl, differentiation_method=diff)
t = np.linspace(0.0, len(Xa)-1, num=len(Xa))
t2 = np.linspace(0.0, len(t)+pred-1, num=len(t)+pred)           
model.fit(Xa, t, u=u[:len(Xa),:])

model.print()
u0 = Xa[0,:]
sim = model.simulate(u0, t2, u=u[:,:])

for i in range(num_variables):
    plt.plot(sim[:, i], label=f'Simulated Variable {i+1}')
    plt.plot(X1[:len(sim), i], '--', label=f'Original Variable {i+1}')

plt.title("SINDy Simulated vs. Original Data")
plt.xlabel("Time Index")
plt.ylabel("Variable Value")
plt.legend()
plt.show()

[Jacob-Stevens-Haas]

There's two issues here. First, your system doesn't really have any controls. Controls are typically variables that aren't related to the variable of interest, or at least aren't well understood. Calculating the controls from the SVD of your variable of interest is not going to be useful.

Secondly, even if your control matrix C1 was not related to Xa, and you wanted to use the SVD, you would choose full_matrices=False. To understand more, I'd have to point you to learn more about the SVD.

Finally, StableLinearSR3 is wonky to begin with. So use at your own risk.

[anfh27]

It's not what's happening with the controls in my real system. Pysindy should (and does) provide a perfect fit for controls = state[1:], so it just makes it easier to see what else is going on.

Regressing on Full_matrices = True gives less variance in the coefficients over time in the real system. Changing to =False does solve the problem here, though unfortunately it doesn't for my real system.
Generally speaking, full_matrices = 1 reliably produces a very good fit on historical data. Perhaps it is overfitting despite the sparse regression? Though, changing nu doesn't seem to have much / any impact.

[anfh27]

Oh, no, I see, thank you. It must be that using full matrices is allowing a kind of false fit on any external shocks in my real system, which is why it's fitting so closely the historical data, but in a way which doesn't allow the best fit of the real dynamics.
Perhaps the fit with full_matrices = 0 is just the best I can achieve with that control.
I will go back to the beginning and check what else might have been incorrect for it to have been producing a similar problem with full_matrices=0.

[Jacob-Stevens-Haas]

Exactly, kind of :)
full_matrices=True is adding what are essentially random vectors to the regression. So it's easier to fit the training data simply because you have more variables to play with. It's purely overfitting.

But again, it's generally expected that you either fit the SVD or the original variables. The SVD is just a linear combination of the variables, so gives you nothing but additional numerical problems in the optimization.

nu typically needs to be played with over a very large range.

[anfh27]

Thanks, that makes sense.
Going back to the beginning, I was able to produce similar errors by delay embedding the state and controls, either by adding too much noise before embedding, or by using the wrong number of delays on state or controls.
I was using cross validation, but median vs mean as the error metric gives very different "optimum" delays, with the median continuing to improve slightly as delays increase, and the mean optimal at low delays but sacrificing the fit. It's hard to know what the "best" metric is for capturing the real relationship between state and controls.

I could only find univariate implementations of mutual information and Cao's FNN in python. Not a pysindy question at all in the end, but if you're aware of any multivariate implementations of other good methods to find tau and delay embedding, please let me know.

[Jacob-Stevens-Haas]

I think I know what you mean by tau, but no, there isn't any good method to discover a classical time delay-differential equation.

For a multivariate delay embedding, there's no problem with emulating the univariate case. If your variables are $x_1(t), x_2(t)$, then a matrix like this should work:

$$\begin{bmatrix} x_1(t_n) & x_2(t_n) & x_1(t_{n-1}) & x_2(t_{n-1}) \\\\\ x_1(t_n) & x_2(t_n) & x_1(t_{n-1}) & x_2(t_{n-1}) \\\\\ \vdots & \vdots & \vdots & \vdots \\\\\ x_1(t_{2}) & x_2(t_2) & x_1(t_1) & x_2(t_1) \end{bmatrix}$$

[anfh27]

That's a shame. Thanks for the help

[Jacob-Stevens-Haas]

you're welcome! And if you find a paper that does something SINDy like for finding tau, I'm happy to countenance adding it in.

[anfh27]

I found that Julia has a few Universal Embedding algorithms which work on multivariate timeseries:
https://github.com/JuliaDynamics/DelayEmbeddings.jl

It should be possible to make the state space then bring it to pysindy to solve, if any work well.
They're not quite what I was thinking of - I had read the HAVOK paper so was imagining multiple successive delays of an unspecified number. But here they pick just a few that accurately reconstructs the attractor - perhaps that can also help with the overfitting.
https://juliadynamics.github.io/DelayEmbeddings.jl/stable/unified/

Pecuzal is in python but apparently much slower
https://pypi.org/project/pecuzal-embedding/

[Jacob-Stevens-Haas]

Thanks, appreciate the links. It'll take a bit of time to read through them, but input embedding, whether via delay embedding, autoencoders, or other ways is something we should look into adding.