input: Write a function to perform gradient descent for linear regression using numpy.
  assert np.allclose(gradient_descent(np.array([[1, 2], [3, 4], [5, 6]]), np.array([1, 2, 3]))[0][-1], np.array([0.091, 0.115]), rtol=1e-2)
  assert np.allclose(gradient_descent(np.array([[1, 1], [1, 2], [2, 2]]), np.array([1, 2, 3]))[0][-1], np.array([0.428, 0.714]), rtol=1e-2)
  assert np.allclose(gradient_descent(np.array([[1, 0], [0, 1], [1, 1]]), np.array([1, 2, 3]))[0][-1], np.array([1, 2]), rtol=1e-2)
output: def gradient_descent(X, y, learning_rate=0.01, iterations=1000):
      import numpy as np
      m, n = X.shape
      theta = np.zeros(n)
      X = np.array(X)
      y = np.array(y)
      theta_history = []
      cost_history = []
      
      for _ in range(iterations):
          predictions = np.dot(X, theta)
          errors = predictions - y
          gradient = np.dot(X.T, errors) / m
          theta = theta - learning_rate * gradient
          
          cost = (1 / (2 * m)) * np.sum(errors ** 2)
          cost_history.append(cost)
          theta_history.append(theta.copy())
      
      return np.array(theta_history), np.array(cost_history)