[GRASS-dev] Compute mahalanobis distance using Scipy

[Glynn Clements]

Paulo van Breugel wrote:

> > The main speed-up will come from "inlining" distance.mahalanobis(),
> > which is essentially:
> >
> >     delta = u - v
> >     m = np.dot(np.dot(delta, VI), delta)
> >     return np.sqrt(m)
> >
> > np.dot(v, m) is equivalent to np.sum(v * m,axis=0), and np.dot(u,v) is
> > equivalent to np.sum(u * v), so the second line above is equivalent to
> >
> >     m = np.sum(np.sum(delta[:,None] * VI * delta[None,:], axis=-1),
> > axis=-1)
> 
> 
> Smart! I can follow the logic, but I am not sure how to solve the problem
> below:
> 
>    Traceback (most recent call last):
>     File "<stdin>", line 1, in <module>
>     ValueError: operands could not be broadcast together with shapes
> (3,1,77,78) (3,3)
> 
> Which refers to the different dimensions of delta and VI?

The first version of the code (which is quoted above) will only work
with 1-D arrays. It's just a fairly direct translation of the existing
distance.mahalanobis() function, given to explain how it can then be
extened to 3-D arrays (i.e. a 2-D array of 1-D vectors).

The second version:

> >     VI = np.linalg.inv(covar)
> >     delta = dat_ref - stat_mean[:,None,None]
> >     m = np.sum(np.sum(delta[:,:,:,None] * delta[:,:,None,:] * VI[None,None,:,:],axis=-1),axis=-1)
> >     stat_mah = np.sqrt(m)

should work with delta being a 3-D array.

-- 
Glynn Clements <glynn at gclements.plus.com>