Linear Regression in Python using gradient descent
That could be due to many different reasons. The most important one is that your cost function might be stuck in local minima. To solve this issue, you can use a different learning rate or change your initialization for the coefficients.
There might be a problem in your code for updating weights or calculating the gradient.
However, I used both methods for a simple linear regression and got the same results as follows:
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.datasets import make_regression


# generate regression dataset
X, y = make_regression(n_samples=100, n_features=1, noise=30)



def cost_MSE(y_true, y_pred):
    '''
    Cost function
    '''
    # Shape of the dataset
    n = y_true.shape[0]
    
    # Error 
    error = y_true - y_pred
    # Cost
    mse = np.dot(error, error) / n
    return mse


def cost_derivative(X, y_true, y_pred):
    '''
    Compute the derivative of the loss function
    '''
    # Shape of the dataset
    n = y_true.shape[0]
    
    # Error 
    error = y_true - y_pred
    
    # Derivative
    der = -2 / n * np.dot(X, error)

    return der


# Lets run an example

X_new = np.concatenate((np.ones(X.shape), X), axis = 1)
learning_rate = 0.1
X_new_T = X_new.T
n_iters = 100
mse = []

#initialize the weight vector
alpha = np.array([0, np.random.rand()])

for _ in range(n_iters):
    
    # Compute the predicted y
    y_pred = np.dot(X_new, alpha)
    
    # Compute the MSE
    mse.append(cost_MSE(y, y_pred))
    
    # Compute the derivative
    der = cost_derivative(X_new_T, y, y_pred)
    
    # Update the weight
    alpha  -= learning_rate * der
alpha

for the gradient descent the coefficients were:
array([-3.36575322, 28.06370831])

Here is the code for closed-form solution:
np.dot(np.linalg.inv(np.dot(X_new_T,X_new)), np.dot(X_new_T, y))

And the coefficients for the closed-form solution:
array([-3.36575322, 28.06370831])

As the coefficients are equal, the RMSE,  MSE, R2 are equal.