title: Multivariate Linear Regression in Python does not find the optimal solution
body: <p>I'm in the progress of implementing my own LinearRegression class from scratch in Python. I am using gradient descent to get the weights and biases for the given X, but even the simplest example does not give a satisfactory output.</p>
  <p>For example:</p>
  <pre class="lang-py prettyprint-override"><code>X = np.array([[1, 1], [1, 2], [2, 2], [2, 3]])
  y = np.dot(X, np.array([1, 2])) + 3 # [6, 8, 9, 11]
  </code></pre>
  <p>gives me coefficients, bias: [1.66112224 2.26236229], 1.3525534054584978, which is clearly not correct. I have tried many different implementations, but I cannot get the correct answer, unfortunately. My R-squared is 0.9238, but I think 1 should be possible, right?</p>
  <p>My class right now:</p>
  <pre class="lang-py prettyprint-override"><code>class LinearRegression:
      def __init__(self):
          self.weights = None
          self.bias = None
  
      def fit(self, X , y, learning_rate=0.001, epochs=1000):        
          # Initialize weights and bias
          if self.weights is None:
              self.weights = np.zeros(X.shape[1])
              self.bias = 0.0
  
          # Perform gradient descent
          for _ in range(epochs):
              y_pred = self.predict(X)
              
              # Calculate the partial derivatives
              partial_w = (1 / X.shape[0]) * (2 * np.dot(X.T, (y_pred - y)))
              partial_b = (1 / X.shape[0]) * (2 * np.sum(y_pred - y))
  
              self.weights -= learning_rate * partial_w
              self.bias -= learning_rate * partial_b
          
          return self
  
      def predict(self, X):
          return X.dot(self.weights) + self.bias
  </code></pre>