Train error vs. Test error in linear regression by samples analysis I have run a multivariate linear regression model on a small set of about 3500 samples. While the model's error is as large as expected, I also ran a bias vs. variance analysis by comparing the train set error vs. the test set error using different sample sizes. I was expecting something like this:

But instead I found that the train error doesn't plateau at any point.

The code that generates this graph is the following:
def loss_by_sample(X, Y, testX, testY, learning_rate):
    samples_list = list()
    train_loss_list = list()
    test_loss_list = list()
    m = X.shape[0]
    for i in range(50, len(X), 50):
        samples_list.append(i)
        weights, loss = gradient_descent(learning_rate, X[:i], Y[:i])

        X2 = tf.concat([tf.ones([X.shape[0], 1]), X], 1)
        train_loss = (1/(2 * m) * tf.tensordot(tf.transpose(h(X2[:i], weights) - Y[:i]), (h(X2[:i], weights) - Y[:i]), axes=1))[0][0]
        train_loss_list.append(train_loss)

        X2 = tf.concat([tf.ones([testX.shape[0], 1]), testX], 1)
        test_loss = (1/(2 * m) * tf.tensordot(tf.transpose(h(X2, weights) - testY), (h(X2, weights) - testY), axes=1))[0][0]
        test_loss_list.append(test_loss)
    // plot train_loss, test_loss

I have tried different learning rates and the one I picked is the one that minimizes the loss using mean squared error:

(values to the right of the x axis after the last datapoint are inf.
Is there something I can conclude from this? Any reason why the train error can surpass the test error?
It occurs because your data is noisy. Suppose that your model is
$$ y = X\beta + noise. $$
Suppose that $\beta$ is recovered exactly. Since the noise is present, our model will NOT predict accurately, it will be always some noise. In other words, think that any time we predict the values, there is some random perturbation. The amount of perturbation is of the same order each time. When we compute the error, we square these perturbations, sum them up, and then take the square root, i.e., compute the norm
In your example, the size of the training set is growing, and, therefore, the test error WITHOUT averaging is growing as well -- you sum up squares of the perturbation. The size of the test set is steady, so the norm of error is about the same there.
Solution: use appropriate averaging to both train and test set.
train_loss = (1/np.sqrt(i) * tf.tensordot(tf.transpose(h(X2[:i], weights) - Y[:i]), (h(X2[:i], weights) - Y[:i]), axes=1))[0][0]
test_loss = (1/(2 * X2.shape[0]) * tf.tensordot(tf.transpose(h(X2, weights) - testY), (h(X2, weights) - testY), axes=1))[0][0]

P.S. X2 in my formulas are as defined in your cell, accordingly. I would suggest choosing a different name for the variable.