Suspected code error in root-function when optimizing with newton_krylov() from scipy

Question

I'm trying to optimize the following equation using the newton_krylov function from scipy:

## This function is what will be used in the newton_krylov optimization process
def residual_function(all_vars):
    """
    Calculates the residual of the weight matrix equation.

    Args:
        all_vars: A combined vector of W elements and lambda_vals.

    Returns:
        A NumPy array representing the element-wise residual of the equation.
    """

    n, m = w_init.shape  # Determine dimensions of W from w_init (this is defined before the   residual_function in my code

    W = all_vars[:n*m].reshape(n, m) ## extracts the first n*m elements from all_vars and assigns to W, reshaping to matrix

    lambda_vals = all_vars[n*m:] ## extracts remaining elements and assigns to lambda_vals

    exp_term = np.exp(-d1 - lambda_vals[0] * Owner - lambda_vals[1] * NH_White) - W ## Owner and NH_white are known column vectors of the same length as W

    residual = d1 * exp_term

    return residual

However, when running:

## set up dual values 
## Combine initial guesses for lambda and gamma
all_lambda = np.array([lambda_init, gamma_init]) ## lambda_init and gamma_init are prespecified as 1 elsewhere in my program

## Combine all necessary inputs, including flattened w_init matrix
all_vars = np.concatenate((w_init.flatten(), all_lambda))

test = residual_function(all_vars = all_vars)

print(test) 
## Solve using newton_krylov with residual function
sol = newton_krylov(residual_function, all_vars)

I'm told the the Jacobian cannot be calculated because the Hessian isn't square.

In this case W is a 23,890 element matrix and I have two duals -- lambda and gamma. My understanding is that the Hessian in this case should be a 23,892 by 23,892 square matrix; however, I the error message states that the matrix is 23,890 x 23,892. I suspect this means that my residual_function() is incorrectly specified.

What I'd expect to be produced in this code is optimal values of lambda and gamma, so a two element array. Any help with this issue would be greatly appreciated.