Problem while increasing the polynomial degree in the polynomial library

[losataschu]

Hello,

I am trying to apply the SINDy algorithm to the following ODE system:
A' = -A
B' = A - 0.5 * B
C' = 0.5 * B
with initial values (A_0, B_0, C_0) = (1,0,0):

# Parameters
dt = 0.005
T = 5
t = np.arange(0, T + dt, dt)
k1 = 1
k2 = 0.5

# reaction model
def rxnrate(t, c, k1=k1, k2=k2):
    A, B, C = c
    dA_dt = -k1 * A
    dB_dt = k1 * A - k2 * B
    dC_dt = k2 * B
    return [dA_dt, dB_dt, dC_dt]

# initial value
x0 = [1.0, 0.0, 0.0]
# numerical result
result = integrate.solve_ivp(
    fun=rxnrate,
    t_span=(t[0], t[-1]),
    y0=x0,
    t_eval=t,
    args=(k1, k2))
x = result.y.T

# noise
np.random.seed(42)  
noise_level = 0.01  
noise = noise_level * np.random.randn(*x.shape)
x_noisy = x + np.random.normal(0,0.01,size=x.shape)

After solving it numerically and adding some noise (0.01 variance), I try to win back the original model using PySINDy and the noisy data. The problem arises when I set the degree of the polynomial library to 2. Even with normalized columns, the resulting ODE system is nowhere close to the original. That doesn't happen if the degree is 1. A second difference between degree 1 and 2 is the conditioning number of the matrix Theta, which also increases a lot in the latter case.

sfd = SmoothedFiniteDifference(smoother_kws={'window_length': 5})
feature_library = ps.PolynomialLibrary(degree=1, include_bias=False)
feature_library2 = ps.PolynomialLibrary(degree=2, include_bias=False)
optimizer = ps.STLSQ(threshold=0.5, normalize_columns=True, max_iter=50)

x_train = x_noisy[0:400].copy()

model_A1s = ps.SINDy(
    differentiation_method=sfd,
    feature_library=feature_library,
    optimizer=optimizer,
)

model_A1s.fit(x_train, t=dt, feature_names=["A","B","C"])
print("Conditioning number poly degree 1:", np.linalg.cond(model_A1s.optimizer.Theta_))
model_A1s.print()

model_A2s = ps.SINDy(
    differentiation_method=sfd,
    feature_library=feature_library2,
    optimizer=optimizer,
)

model_A2s.fit(x_train, t=dt, feature_names=["A","B","C"])
Theta_A2s = model_A2s.optimizer.Theta_
print("Conditioning number poly degree 2:", np.linalg.cond(Theta_A2s))
model_A2s.print()

The coefficients look like this:

Conditioning number poly degree 1: 4.478519350721915
(A)' = -0.915 A + -0.104 B +  0.068 C
(B)' =  1.054 A + -0.524 B
(C)' =  0.056 A +  0.415 B +  0.067 C

Conditioning number poly degree 2: 801.7100259304184
(A)' = -0.584 A + -0.269 A^2 + -0.753 A B
(B)' = -0.129 A +  0.071 B +  1.310 A^2 + -0.676 B C
(C)' =  0.653 B +  1.996 A C + -0.836 B^2

Would you have any idea about the cause behind this behavior? Does the ODE system with initial conditions need to fulfill some conditions in order to apply PySINDy successfully? How would you tackle it?

[losataschu]

Update:

I replicated the example with the Lorenz system https://pysindy.readthedocs.io/en/stable/examples/5_differentiation/example.html in my laptop and I was able to get some results (the code shown in the website is a bit outdated, but it was easy to adapt it).

The conditioning number of the matrices Theta in both cases, clean and noisy data, were around 1000, so it appears to me that the problem must lie somewhere else, but I'm not able no identify where yet.

I'd be grateful for your help or guidance in this regard. Thank you!

[Jacob-Stevens-Haas]

There's a few things that seemed to be going on when I ran your code.

1. It looks like sfd.smoothed_x_ is still pretty noisy. Increasing the window to 100 makes it look much better and improves both models performance.
2. Since you're adding a lot of terms to the feature library, you might need a stronger regularizer (alpha=1.0 or more). This didn't fully fix the problem
3. Condition number matters a lot in sparse regression, because it's not just about numerically distinguishing the features. The guarantees in the original L-1 regularizer paper only applied with low condition numbers (or more correctly, in the presence of a Null Space Property and Restriced Isometry Property). The effect varies between optimizers, however. MIOSR does much better here, e.g.

[losataschu]

Thank you very much for the help. I'll check on MIOSR and its implementation.

[losataschu]

Hello Jacob,

thank you again for your suggestions. The optimizer MIOSR worked quite well, both target and group sparsity identified the relevant terms. It still seems a bit forced that I first need to determine how many coefficients will be non zero.
When you ran this code, do you remember if the result was close enough to the original model?

Best,
Lorenzo

[Jacob-Stevens-Haas]

I don't remember. However, certain datasets actually work better with some noise. A good example is an object in free fall:

$$ \dot x = v; ~~~~~ \dot v = -a x$$

A polynomial library in x and v has condition number infinity: a linear combination of terms (x and v**2) will give you the total energy. Since total energy is conserved, it is constant, which is collinear with a bias term. That system is also linear, and the problem only arises when the polynomial library grows to degree 2.

I don't have a bulletproof way to detect that a system lies on a linear manifold of the feature library, other than (a) verify on noiseless data at high resolution to rule out differentiation, (b) try with smaller or larger feature libraries, and (c) check the condition number. If the error grows when noise goes away and if the problem only appears when the feature library adds some terms, the problem may be that solutions lie on a manifold.