How do I fit a multi-linear regression model with a Loss Function using the scipy.minimize function?

Question

I'm currently learning multilinear regression and I'm stuck on this question that wants me to define a function that will "Minimize the average loss calculated from using different theta vectors, and estimate the optimal theta for the model".

from scipy.optimize import minimize
    
    def l1(y, y_hat):
        return np.abs(y - y_hat)
    
    def l2(y, y_hat):
        return (y - y_hat)**2
    
    def minimize_average_loss(loss_function, model, X, y):
        """
        Minimize the average loss calculated from using different theta vectors, and 
        estimate the optimal theta for the model.
        
        Parameters
        -----------
        loss_function: either the squared or absolute loss functions defined above
        model: the model (as defined in Question 1b)
        X: a 2D dataframe of numeric features (one-hot encoded)
        y: a 1D vector of tip amounts
        
        Returns
        -----------
        The estimate for the optimal theta vector that minimizes our loss
        """
        
        ## Notes on the following function call which you need to finish:
        # 
        # 0. The first '...' should be replaced with the average loss evaluated on 
        #       the data X, y using the model and appropriate loss function.
        # 1. x0 are the initial values for THETA.  Yes, this is confusing
        #       but optimization people like x to be the thing they are 
        #       optimizing. Replace the second '...' with an initial value for theta,
        #       and remember that theta is now a vector. DO NOT hard-code the length of x0;
        #       it should depend on the number of features in X.
        # 2. Your answer will be very similar to your answer to question 2 from lab 7.
        ...
        return minimize(lambda theta: ..., x0=...)['x']
        # Notice above that we extract the 'x' entry in the dictionary returned by `minimize`. 
        # This entry corresponds to the optimal theta estimated by the function.
    
    minimize_average_loss(l2, linear_model, one_hot_X, tips)

For context, my linear model is defined as this:

def linear_model(thetas, X):
    """
    Return the linear combination of thetas and features as defined above.
    
    Parameters
    -----------
    thetas: a 1D vector representing the parameters of our model ([theta1, theta2, ...])
    X: a 2D dataframe of numeric features
    
    Returns
    -----------
    A 1D vector representing the linear combination of thetas and features as defined above.
    """
    return np.dot(X, thetas)

Currently, I have:

def minimize_average_loss(loss_function, model, X, y):
    return minimize(lambda theta: loss_function(y, linear_model(theta, X)), x0= [0.0, 0.0])['x']

Does anyone know how I'm supposed to do this? Thanks!