Effectively solve an overdetermined nonlinear equation system using fitted data in python-CodePudding

I think it’s the easiest way to describe my problem with a small example. I have this data which is my input data. I have 3 LEDs and each LED is represented by 4 color-parts (value 1 to 4 in each array). If I increase the intensity of the LEDs (in this example from 10% to 30%) the color-parts change in different ways.

LED1_10 = np.array([1.5, 1, 0.5, 0.5])
LED1_20 = np.array([2.5, 1.75, 1.2, 1.2])
LED1_30 = np.array([3, 2.3, 1.7, 1.7])

LED2_10 = np.array([0.2, 0.8, 0.4, 0.4])
LED2_20 = np.array([0.6, 1.6, 0.5, 0.5])
LED2_30 = np.array([1.0, 2.0, 0.55, 0.55])

LED3_10 = np.array([1, 0.1, 0.4, 0.4])
LED3_20 = np.array([2.5, 0.8, 0.9, 0.9])
LED3_30 = np.array([3.25, 1, 1.3, 1.3])

The column elements of the arrays belong together. So if I bring LED1 from 10% to 30% the value in column 1 rises from 1.5 to 2.5 and then to 3. I want to find a polynomial for the rise of each of the LED color-parts, so I rearrange the data and use a polynomial fit to get the equations which describe the way the values are rising for each LED.

### Rearrange the values
LED1 = np.stack((LED1_10, LED1_20, LED1_30)).T
LED2 = np.stack((LED2_10, LED2_20, LED2_30)).T
LED3 = np.stack((LED3_10, LED3_20, LED3_30)).T

### Create x-vectro
x = np.array([10,20,30])

### Polynomal fits
Fit_LED1 = []
for i in range(len(LED1)):
    z = np.polyfit(x, LED1[i], 2)
    Fit_LED1.append(z)

Fit_LED2 = []
for i in range(len(LED2)):
    z = np.polyfit(x, LED2[i], 2)
    Fit_LED2.append(z)
    
Fit_LED3 = []
for i in range(len(LED3)):
    z = np.polyfit(x, LED3[i], 2)
    Fit_LED3.append(z)

Now I want to generate light of a specific color mixing together the light of each of the 3 different LEDs. Therefor I need to find out which intensity I need to use from each of the LEDs to get the best possible result. Color-parts 1-4 are represented by the solution vector: b = [7, 8, 2, 5] I do this solving the overdetermined nonlinear equation system like this:

def f(x):
    x1, x2, x3 = x
    return np.asarray(((-2.50000000e-03*x1**2   1.75000000e-01*x1   -5.91091254e-15)   (-2.03207837e-18*x2**2   4.00000000e-02*x2   -2.00000000e-01)   (-0.00375*x3**2   0.2625*x3   -1.25),
                       (-0.001*x1**2   0.105*x1    0.05)   (-0.002*x2**2    0.14*x2   -0.4)   (-0.0025*x3**2    0.145*x3   -1.1),
                       (-0.001*x1**2   0.1*x1   -0.4 )   (-0.00025*x2**2    0.0175*x2    0.25)   (-0.0005*x3**2    0.065*x3   -0.2), 
                       (-0.001*x1**2   0.1*x1   -0.4 )   (-0.00025*x2**2    0.0175*x2    0.25)   (-0.0005*x3**2    0.065*x3   -0.2)))

def system(x,b):
    return (f(x)-b)

b = [7, 8, 2, 5]
x = scipy.optimize.leastsq(system, np.asarray((1,1,1)), args=b)[0]

I used the polynomials I got from the fit for each LED and added them to each other to get a equation for each of the four color parts. This would give the same result and maybe is a bit easier to read:

def g(x):
    x1, x2, x3 = x
    return np.asarray(((Fit_LED1[0][0]*x1**2   Fit_LED1[0][1]*x1   Fit_LED1[0][2])   (Fit_LED2[0][0]*x1**2   Fit_LED2[0][1]*x1   Fit_LED2[0][2])   (Fit_LED3[0][0]*x1**2   Fit_LED3[0][1]*x1   Fit_LED3[0][2]),
                       (Fit_LED1[1][0]*x1**2   Fit_LED1[1][1]*x1   Fit_LED1[1][2])   (Fit_LED2[1][0]*x1**2   Fit_LED2[1][1]*x1   Fit_LED2[1][2])   (Fit_LED3[1][0]*x1**2   Fit_LED3[1][1]*x1   Fit_LED3[1][2]),
                       (Fit_LED1[2][0]*x1**2   Fit_LED1[2][1]*x1   Fit_LED1[2][2])   (Fit_LED2[2][0]*x1**2   Fit_LED2[2][1]*x1   Fit_LED2[2][2])   (Fit_LED3[2][0]*x1**2   Fit_LED3[2][1]*x1   Fit_LED3[2][2]),
                       (Fit_LED1[3][0]*x1**2   Fit_LED1[3][1]*x1   Fit_LED1[3][2])   (Fit_LED2[3][0]*x1**2   Fit_LED2[3][1]*x1   Fit_LED2[3][2])   (Fit_LED3[3][0]*x1**2   Fit_LED3[3][1]*x1   Fit_LED3[3][2])))
def system(x,b):
    return (f(x)-b)

b = [5, 8, 4, 12]
x = scipy.optimize.leastsq(system, np.asarray((1,1,1)), args=b)[0]

Now my problem is that I need to type each of the functions separately which is a lot of work, especially because my real-world application consists of 40 LEDs and more than 1000 color-parts for each of the LEDs. Is there an easier and more efficient way to define the equations for the equations system rather than typing each of the separately like I did here?

def g(x):
    x1, x2, x3 = x
    return np.asarray(((Fit_LED1[0][0]*x1**2   Fit_LED1[0][1]*x1   Fit_LED1[0][2])   (Fit_LED2[0][0]*x1**2   Fit_LED2[0][1]*x1   Fit_LED2[0][2])   (Fit_LED3[0][0]*x1**2   Fit_LED3[0][1]*x1   Fit_LED3[0][2]),
                       (Fit_LED1[1][0]*x1**2   Fit_LED1[1][1]*x1   Fit_LED1[1][2])   (Fit_LED2[1][0]*x1**2   Fit_LED2[1][1]*x1   Fit_LED2[1][2])   (Fit_LED3[1][0]*x1**2   Fit_LED3[1][1]*x1   Fit_LED3[1][2]),
                       (Fit_LED1[2][0]*x1**2   Fit_LED1[2][1]*x1   Fit_LED1[2][2])   (Fit_LED2[2][0]*x1**2   Fit_LED2[2][1]*x1   Fit_LED2[2][2])   (Fit_LED3[2][0]*x1**2   Fit_LED3[2][1]*x1   Fit_LED3[2][2]),
                       (Fit_LED1[3][0]*x1**2   Fit_LED1[3][1]*x1   Fit_LED1[3][2])   (Fit_LED2[3][0]*x1**2   Fit_LED2[3][1]*x1   Fit_LED2[3][2])   (Fit_LED3[3][0]*x1**2   Fit_LED3[3][1]*x1   Fit_LED3[3][2])))

I hope I was able to make my problem clear and would be very thankful if anyone could help me solving this task.

Thank you very much in advance :)

CodePudding user response：

There is a pattern in your equations which you can vectorise. First of all, gather all the LED fits in a 3D array.

fits = np.array([Fit_LED1, Fit_LED2, Fit_LED3])

And then define g(x) as

def g(x):
    X = np.array([x**2, x, np.ones_like(x)]).T
    return np.sum(fits * X[:,None], axis=(0, 2))

You can also confirm the result is correct with np.isclose(f(x), g(x)).

Of course you should do the same to LED1, LED2, etc so you don't have to hardcode Fit_LED1, etc. Just put everything in a 3d array and loop for each LED index.

LEDs = np.array([LED1, LED2, LED3])
fits = [
    [np.polyfit(np.array([10, 20, 30]), LEDs[i,j], 2) for j in range(LEDs.shape[1])]
    for i in range(LEDs.shape[0])
]
fits = np.array(fits)