Numba GuFunc giving incorrect output

Question

I have a python function that follows:

    @njit(cache=True, fastmath=True, parallel=True)
    def min_var_opt(covariance_matrix: np.ndarray) -> np.array:
        n = covariance_matrix.shape[0]

        # Compute the inverse of the covariance matrix
        inverse_covariance = np.linalg.inv(covariance_matrix)
    
        ones_vector = np.ones((n, 1))
    
        denominator = np.dot(np.dot(ones_vector.T, inverse_covariance), ones_vector)
    
        weights = np.dot(inverse_covariance, ones_vector) / denominator
        
        return weights
    
    covariance_matrix = np.array(
            [
                [77.25607201, 4.95549419, -2.1582784],
                [4.95549419, 73.87388998, -3.76609601],
                [-2.1582784, -3.76609601, 259.46734795],
            ]
        )
    
    weights = min_var_opt(covariance_matrix).T
    print(f"compiled output:\t{weights}")

This produces the correct output: compiled output: [[0.41740475 0.4411879 0.14140735]]

Now, I would like to generalize this over a series of covariance matrices:

    @guvectorize([(float64[:, :], float64[:])], "(n, n) -> (n)", nopython=True, cache=True)
    def minimum_variance_optimization(covariance_matrix: np.ndarray, weights: np.array):
        n = covariance_matrix.shape[0]
    
        # Compute the inverse of the covariance matrix
        inverse_covariance = np.linalg.inv(covariance_matrix)
    
        ones_vector = np.ones((n, 1))
    
        denominator = np.dot(np.dot(ones_vector.T, inverse_covariance), ones_vector)
    
        weights = np.dot(inverse_covariance, ones_vector) / denominator
    
    cov = np.array([covariance_matrix])
    vectorized_weights = minimum_variance_optimization(cov)
    print(f"vectorized output once:\t{vectorized_weights}")

This gives the identical output: vectorized output once: [[0.41740475 0.4411879 0.14140735]]

This is the expected output. The issue seems to be when i apply it with more than once inner matrix.

When I try:

    cov = np.tile(covariance_matrix, (5, 1, 1))
    vectorized_weights = minimum_variance_optimization(cov)
    print(f"vectorized output repeated:\t{vectorized_weights}")

I get the following incorrect output:

    vectorized output repeated:[[4.66464261e-310 0.00000000e+000 0.00000000e+000]
    [4.66394411e-310 3.95252517e-323 0.00000000e+000]
    [3.14015482e+179 2.36897524e+261 2.93013167e-002]
    [7.10264216e+083 1.32448107e+007 3.12881791e+011]
    [6.73475158e+107 4.17959458e+025 5.71862049e-310]]

According to this answer, I need to change the assignment above to be weights[:] = np.dot(inverse_covariance, ones_vector) / denominator to override the input array, but when I try this, numba fails to compile, and I don't understand how or why but it is writing a result with the original code above

I expect the output to be the same as the first vector repeated 5 times. I must be invoking the the gufunc incorrectly, but I cannot find any examples that match this use case. Can someone help me understand what is happening?

maciek97x · Accepted Answer · 2023-07-12T07:17:48.967

You need to use for loop for the first dimension of input array. Here is fixed implementation:

@nb.guvectorize([(nb.float64[:, :, :], nb.float64[:, :])], "(k, n, n) -> (k, n)", nopython=True, cache=True)
def minimum_variance_optimization(covariance_matrix: np.ndarray, weights: np.ndarray):
    for k in range(covariance_matrix.shape[0]):
        n = covariance_matrix.shape[1]

        # Compute the inverse of the covariance matrix
        inverse_covariance = np.linalg.inv(covariance_matrix[k, ])

        ones_vector = np.ones((n, 1))

        denominator = np.dot(np.dot(ones_vector.T, inverse_covariance), ones_vector)

        weights[k] = (np.dot(inverse_covariance, ones_vector) / denominator)[:, 0]

Results:

cov = covariance_matrix[np.newaxis, ]
vectorized_weights = minimum_variance_optimization(cov)
print(f"vectorized output once:\t{vectorized_weights}")

vectorized output once: [[0.41740475 0.4411879  0.14140735]]

cov = np.tile(covariance_matrix, (5, 1, 1))
vectorized_weights = minimum_variance_optimization(cov)
print(f"vectorized output repeated:\t{vectorized_weights}")

vectorized output repeated: [[0.41740475 0.4411879  0.14140735]
 [0.41740475 0.4411879  0.14140735]
 [0.41740475 0.4411879  0.14140735]
 [0.41740475 0.4411879  0.14140735]
 [0.41740475 0.4411879  0.14140735]]

Thanks for replying, your solution worked! I guess I never saw any example of higher dimension arrays, but I'll keep in mind that I need to loop over the entire output/input. You saved me a lot of time, thanks again! — user3299166, Jul 12 '23 at 16:44

Numba GuFunc giving incorrect output

1 Answers1