Руководство по более эффективному применению метода наименьших квадратов

Я реализовал программу, которая выполняет обычные методы наименьших квадратов в необработанном питоне (но с использованием массивов numpy, а не списков). Я бы хотел, чтобы мой код был раскритикован сообществом.

Что является более быстрым / более эффективным способом реализации обычного метода наименьших квадратов, чем то, что я сделал в моем фрагменте кода ниже?

Предполагается, что существует одно измерение функции. (AT * A равно 2 X 2)

Моя цель — реализовать это на C, а затем распараллеливать части в CUDA C

mat = np.array([
    [1.,1.],
    [1.,2.],
    [1.,3.],
])

arr = np.array([
                [1., 2., 4.]
]).T


def least_squares(mat, arr):
    #matrix for A.T * A
    new_arr = np.array([
                [0.,0.],
                [0.,0.]
    ])
    #Right hand side vector
    new_right = np.array([
            [0., 0.]
    ]).T

    #Solution Vector
    x_hat = np.array([
            [0., 0.]
    ]).T
    #Calculate A.T * A
    for i in range(3):
        new_right[0][0] += mat[i][0] * arr[i]
        new_right[1][0] += mat[i][1] * arr[i]
        for j in range(2):
            new_arr[0][j] += mat[i][0] * mat[i][j]
            new_arr[1][j] += mat[i][1] * mat[i][j]

    #Perform Elimination on A.T * A & right side b        
    row_multiple = new_arr[1][0]/new_arr[0][0]
    new_arr[1] -= (row_multiple)*new_arr[0]
    new_right[1,0] -= row_multiple*new_right[0,0]

    
    #Solve for x-hat
    x_hat[1,0] = new_right[1,0] / new_arr[1,1]
    x_hat[0,0] = (new_right[0,0] - new_arr[0,1]*(new_right[1,0] / new_arr[1,1])) / new_arr[0,0]
    return x_hat