Squared Mahalanobis distance function in Python returning array - why?

Question

The code is:

import numpy as np
def Mahalanobis(x, covariance_matrix, mean):
    x = np.array(x)
    mean = np.array(mean)
    covariance_matrix = np.array(covariance_matrix)
    return (x-mean)*np.linalg.inv(covariance_matrix)*(x.transpose()-mean.transpose())

#variables x and mean are 1xd arrays; covariance_matrix is a dxd matrix
#the 1xd array passed to x should be multiplied by the (inverted) dxd array
#that was passed into the second argument
#the resulting 1xd matrix is to be multiplied by a dx1 matrix, the transpose of 
#[x-mean], which should result in a 1x1 array (a number)

But for some reason I get a matrix for my output when I enter the parameters

Mahalanobis([2,5], [[.5,0],[0,2]], [3,6])

output:

out[]: array([[ 2. ,  0. ],
              [ 0. ,  0.5]])

It seems my function is just giving me the inverse of the 2x2 matrix that I input in the 2nd argument.

Answer 1

You've made the classic mistake of assuming that the * operator is doing matrix multiplication. This is not true in Python/numpy (see http://www.scipy-lectures.org/intro/numpy/operations.html and https://docs.scipy.org/doc/numpy-dev/user/numpy-for-matlab-users.html). I broke it down into intermediate steps and used the dot function

import numpy as np
def Mahalanobis(x, covariance_matrix, mean):

    x = np.array(x)
    mean = np.array(mean)
    covariance_matrix = np.array(covariance_matrix)

    t1 = (x-mean)
    print(f'Term 1 {t1}')

    icov = np.linalg.inv(covariance_matrix)
    print(f'Inverse covariance {icov}')

    t2 = (x.transpose()-mean.transpose())
    print(f'Term 2 {t2}')

    mahal = t1.dot(icov.dot(t2))

    #return (x-mean)*np.linalg.inv(covariance_matrix).dot(x.transpose()-mean.transpose())
    return mahal

#variables x and mean are 1xd arrays; covariance_matrix is a dxd matrix
#the 1xd array passed to x should be multiplied by the (inverted) dxd array
#that was passed into the second argument
#the resulting 1xd matrix is to be multiplied by a dx1 matrix, the transpose of 
#[x-mean], which should result in a 1x1 array (a number)


Mahalanobis([2,5], [[.5,0],[0,2]], [3,6])

produces

Term 1 [-1 -1]
Inverse covariance [[2.  0. ]
 [0.  0.5]]
Term 2 [-1 -1]
Out[9]:    2.5