I have a set of N samples (N~10000 to 100000) :
(y_i, P_i)
They sample an unknown function :
y = f(P)
In my case, P is a set of n_p parameters with n_p typically around 10.
My aim is to use the samples to find the polynomial of order 2 which approximate the best my unknown function. The typical result of this would be the (n_p+1)+(n_p+1)*n_p/2 coefficients of the best-fit polynomial.
However, I would like to obtain the solution in the following form :
f(P) = (P-mu)*(C^-1)*(P-mu)^T + K
with mu a vector of dimension n_p, C a (n_p X n_p) symmetric matrix, and K a constant. These mu, C and K are what I'd like to obtain.
Is there a standard implementation somewhere in the Python ecosystem and/or an efficent way to do this ?
Thanks in advance !
Edit : Here is a typical setup, where I create fake samples (P and Y) and by only using them, I would like to recover the original mu, C and K :
import numpy as np
n_p = 3
N = 15
# Generate the "mu" we'll try to recover
mu = np.random.rand(n_p)
# Generate the "C" we'll try to recover
tmp = np.random.rand(n_p,n_p)
C = np.dot(tmp, tmp.T)
# Generate the "K" we'll try to recover
K = np.random.rand()
# Generate the samples
P = np.random.rand(n_p,N)
Y = np.array([np.dot(np.dot((P[:,i]-mu),np.linalg.inv(C)),(P[:,i]-mu))+K for i in range(N)])
Is something like this what you're looking for? I'm not 100% sure it's accurate for your use case (I don't know what additional constraints you can impose when you're fitting to an order 2 polynomial), but it should try to find a U, A (what you called C), and K that minimize the least squares error.
There's probably a more mathematically sound approach to this, but hopefully this at least helps.
Edit: Just noticed I messed up the sign of k. I also needed to make the number of parameters higher. This seems to work super well though. I'm pretty much exactly recovering the original noise.
Also adding sample output.