MLP parameters vs train set size for function fitting: right approach?

Question

(positions, velocities, accelerations) = ($q$, $\dot{q}$, $\ddot{q}$) are generalized coordinates of a Lagrangian system

I'm trying to learn the Lagrangian L($q$, $\dot{q}$) by comparing analytical accelerations and predicted accelerations. These are found by manipulating Euler-Lagrange equations into $\ddot{q}$=f(L)=$(\nabla^{\space}_{\dot{q}}\nabla^T_{\dot{q}}L)^{-1}(\nabla_qL-(\nabla^{\space}_q\nabla_{\dot{q}}^TL)\dot{q})$
In synthesis, this is my forward method:
MLP Input ($q$, $\dot{q}$) (4 scalars, because we consider 2D systems)
MLP Output L($q$, $\dot{q}$) (1 scalar)
Calculate $\ddot{q}$= f(L)

These accelerations are then compared with analytical accelerations given by the same q, qt with an MSE loss function.
This happens pointwise in time, in the sense that the set ($q$, $\dot{q}$, $\ddot{q}$) do not carry time information, they're a triplet (sextuplet, because each has D=2) of scalars.
The predicted accelerations are not that bad, but when calculating predicted $q$, $\dot{q}$ by integration, even a small error propagates to a noticeably wrong trajectory.
Before enduring the painful road of hyperparameter optimization, I want to make sure my intuition on the issue is right, so I can confidently assume my MLP structure is not the main bottleneck. Please correct me on anything I'm missing.

Because I'm trying to fit a scalar function, the MLP can be thought of as an equation in P parameters, where P is the number of trainable params of the net.
Training it and comparing predictions with N targets (analytical $\ddot{q}$, so how many doublets $q$, $\dot{q}$ i'm feeding to the MLP) imposes N restraints, so I can think of it as solving a system of N equations in P parameters.
As such, my intuition tells me I need to have AT LEAST N=P points ($q$, $\dot{q}$) to correctly solve it.
In the following example, I have P = 3361. Feeding it N = 4096 points will be enough to train them. Training with N = 16384 points should be better on one hand (because a larger training set can better generalise the predictions, given some initial conditions), but following the N eq in P variables reasoning, it feels like overdetermining the system. I have tried various configurations of this, and having N >> P seems better at fitting $\ddot{q}$, though never providing correct trajectiories because of the aformentioned propagation of errors in integration.

class LNN(nn.Module):
    def __init__(self):
        super(LNN, self).__init__()
        self.fc1 = nn.Linear(4, 32)
        self.fc2 = nn.Linear(32, 32)
        self.fc3 = nn.Linear(32, 32)
        self.fc4 = nn.Linear(32, 32)
        self.fc5 = nn.Linear(32, 1)
        self.sp = nn.Softplus(beta=0.5)   #beta = 1 is steeper, so derivative = 1 is approached faster 

    def lagrangian(self, x): 
        x = self.sp(self.fc1(x))
        x = self.sp(self.fc2(x))
        x = self.sp(self.fc3(x))
        x = self.sp(self.fc4(x))
        x = self.fc5(x)
        return x     
    
    def forward(self, x):  #   q, qt
        n = x.shape[1]//2  #    2   
        xv = torch.autograd.Variable(x, requires_grad=True)                  
        xv_tup = tuple([xi for xi in x]) 
        tqt = xv[:, n:]    #   qt  

        jacpar = partial(jacobian,  self.lagrangian, create_graph=True)
        hesspar = partial(hessian,  self.lagrangian, create_graph=True)
      
        A = tuple(map(hesspar, xv_tup))    #   
        B = tuple(map(jacpar, xv_tup))            ####################      
                                                  ####THIS IS f(L)####
        multi = lambda Ai, Bi, tqti, n:  torch.pinverse(Ai[n:, n:]) @ (Bi[:n, 0] - Ai[n:, :n] @ tqti) 
        multi_par = partial(multi, n=n)
        tqtt_tup = tuple(map(multi_par, A, B, tqt))           
        tqtt = torch.cat([tqtti[None] for tqtti in tqtt_tup]) 
                                                              
        xt = torch.cat([tqt, tqtt], axis=1)
        xt.retain_grad()
        return xt                        #qt (same as input), qtt 

    def t_forward(self, t, x):
        return self.forward(x)

For the training part, I'm feeding the net analytical trajectories where points are shuffled, with the following parallel in mind:
training on batches of a continuous part of a curve, then going through the next part and so on, is like training an image recognition net by feeding ALL cats first, then ALL dogs etc, slowing down the learning of different features by focusing on the ones of a specific class. Again, this is (i think) true because of the pointwise nature of the training.
Here x = ($q$, $\dot{q}$), xt = ($\dot{q}$, $\ddot{q}$) (it's easier to pair them like this for how the code is built)

x = torch.tensor(anal_solve_ode(x0[:2], x0[2:], t_train)) #trajectories from analytical qtt and init. cond.
xt = torch.tensor(get_xt_anal(x, t_train)).float()  #get analytical accelerations for each q, qt i x
ind = torch.randperm(x.shape[0])
x = x[ind, :]   #shuffle
xt = xt[ind, :]

This method, together with training with trajectories coming from different initial conditions, seems to result in better generalisation.

Any input or experience on the matter?

Answer 1

One point to consider is that main modern neural networks are substantially over-parameterized by "traditional" standards. And that's not necessarily a bad thing and even happens in more traditional settings e.g. with splines as illustrated in this nice discussion of the double-descent phenomenon. So, I would not assume the number of trainable parameters needs to be less than the number of training samples, although that may be true. Additionally, there's a lot of recent work on this kind of set-up (e.g. here), so it's worth researching that.