Newton-Gauss method on python for approximation

Question

I am trying to implement the Newton-Gauss method on a python. Here is my full code:

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

def gauss_newton(X, Y, max_iter=1000, eps=1e-6): 
    P0 = [1, 1, 1]
    J = np.zeros([len(X), len(P0)])

    for i in range(max_iter):
        j1 = 1
        j2 = P0[0]
        j3 = P0[2]*X
        J[:,0] = j1
        J[:,1] = j2
        J[:,2] = j3

        r = Y - (P0[0] + P0[1]*X + P0[2]*X**2)
        t1 = np.linalg.inv(np.dot(J.T, J))
        t2 = np.dot(t1, J.T)
        t3 = np.dot(t2, r)

        P1 = P0 - t3
        t4 = abs(P1-P0)
        if max(t4) <= eps:
            break
        P0 = P1

    return P0[0], P0[1], P0[2]

X = np.asarray([1, 2, 3, 4, 5, 6,])
Y = np.asarray([5, 7, 9, 11, 13, 11])
C1, C2, C3 = gauss_newton(X, Y)
pred = C1*X/(C2+X)
plt.figure(1, figsize=(6,4), dpi=120)
plt.scatter(x=X, y=Y, c='green', marker='o', label='Data')
plt.plot(X, pred, '--m', label='Model')
plt.legend()
plt.show()

Unfortunately, there is a "Singilar matrix" error in the line:

t1 = np.linalg.inv(np.dot(J.T, J))

I completely repeated all the points from the video https://www.youtube.com/watch?v=weuKzGFVQV0&t=109s except the used model - [Y = C1 + C2X + C3X**2], not [Y = C1*X/C2+X]. So the derivatives j1,j2, and also j3 (since there are 3 parameters in this model) are calculated differently. The same is true for the parameter r. Here's how they were calculated for the model in the video:

 j1 = -(X/(P0[1]+X))
 j2 = ((P0[0]*X)/(P0[1]+X)**2)
 r = Y - ((P0[0]*X)/(P0[1]+X))
 # also matrix P0 is initialized as P0 = [1,1]

Can you please explain what might be the problem with my code above?

Answer 1

Not a python problem but a derivatives problem.

Having Y = P0[0] + P0[1]*X + P0[2]*X**2, J is found as partial derivatives with respect to the constants, not the variables.

In this case:

j1 = d(Y)/d(P[0]) = 1
j2 = d(Y)/d(P[1]) = X
j3 = d(Y)/d(P[2]) = X**2

The singular matrix is due to having two rows with the same values (j1=1, j2=1) as both vectors are linearly dependant.