title: Fitting polynomials of different degrees to function with custom loss function
body: <p>I am trying to fit polynomial functions of different degrees to some data I generated and I am doing gradient descent without a library. I am also using a custom loss function that I manually calculated the gradient of (hopefully it is correct). The problem is I am getting infinity and nan for the mean of the squared residuals and I have no idea what I am doing wrong. Please help me.</p>
  
  <pre><code>import numpy as np
  import matplotlib.pyplot as plt
  from sklearn.model_selection import train_test_split
  
  #y(x) = sin(x) + noise
  def F(x_in_rad):
      return np.sin(x_in_rad) + noise*np.random.normal(0,1,N)
  
  noise = 0.5
  N = 50
  #X - N datapoints in radians
  X = np.deg2rad(np.random.normal(0,1,N) * 359)
  Y = F(X)
  
  
  X = np.atleast_2d(X).T
  Y = np.atleast_2d(Y).T
  #split data
  X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=0.33, random_state=42)
  
  #mean of the squared residuals
  def sq_res_mean(Y_pred,Y_real):
      return np.mean((Y_pred - Y_real)**2)
  
  #create design matrix
  def designmatpoly(X,degree):
      X = X[:,0]
      eye = np.ones(X.shape)
      rows = []
      rows.append(eye)
      for i in range(1,degree+1):
          rows.append(X**i)
      return np.stack(rows).T
  
  
  #L2 norm squared; gradient = 2w
  def C(w):
      return np.sum(w**2)
  
  
  def gradientdescent(Amat, y, rate, numiter, lam, deg):
      n, p = Amat.shape
      whistory = []
      w_analytical = np.dot((np.dot(Amat.T,Amat) + lam*np.eye(deg+1, dtype=int))**(-1),np.dot(Amat.T,Y_train))
      losshistory = [] 
      #random weights initialized
      w = np.atleast_2d(np.random.randn(deg+1)).T
      for i in range(numiter): 
          loss = np.square(y - w[0] - np.dot(Amat, w)) + lam*C(w_analytical)
          whistory.append(w)
          losshistory.append(loss)
          grad = np.dot(-2*Amat.T, y - w[0] - Amat.dot(w)) + lam*2*w_analytical
          w = w - rate*grad
      return w, np.asarray(whistory), np.asarray(losshistory)
  
  def model(degree,rate=0.0001, num_iters = 50, lam = 0.5):
      A_test =  designmatpoly(X_test,degree)
      A_train = designmatpoly(X_train,degree)
      wfin, whist, meanlosstrace = gradientdescent(A_train, Y_train, rate, num_iters, lam, degree)
      return wfin, A_test
  
  degrees = []
  sq_res_means = []
  for i in range(1,10):
      wfin, A_test = model(degree=i)
      degrees.append(i)
      Y_pred = np.dot(A_test,wfin)
      sqrm = sq_res_mean(Y_pred,Y_test)
      sq_res_means.append(sqrm)
      print("deg",i,"sq_res_mean",sqrm)
  </code></pre>
  
  <p>Loss function:
  <a href="https://i.stack.imgur.com/e1ayj.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/e1ayj.png" alt="enter image description here"></a></p>