马哈拉诺比斯距离反转协方差矩阵

问题描述 投票:0回答:4

我正在编写一个函数来计算两个向量之间的马氏距离。据我所知,这是使用方程 a'*C^-1*b 实现的,其中 a 和 b 是向量,C 是协方差矩阵。我的问题是,是否有一种有效的方法可以在不使用高斯乔丹消除的情况下找到矩阵的逆,或者没有办法解决这个问题?我正在寻找一种自己完成此操作的方法,而不是使用任何预定义的函数。

我知道 C 是 Hermitian 正定矩阵,那么有什么方法可以在算法上利用这一事实吗?或者是否有一些巧妙的方法来计算马氏距离而不计算协方差的倒数?任何帮助将不胜感激。

***编辑:上面的马氏距离方程是不正确的。它应该是 x'*C^-1*x,其中 x = (b-a),b 和 a 是我们试图找到其距离的两个向量(感谢 LRPurser)。因此,所选答案中提出的解决方案如下:

d=x'*b,其中 b = C^-1*x C*b = x,因此使用 LU 分解或 LDL' 分解求解 b。

c algorithm matrix linear-algebra
4个回答
11
投票

您可以(而且应该!)使用 LU 分解,而不是显式反转矩阵:使用分解求解

C x = b
比计算
C^-1
和乘以向量
b
具有更好的数值属性。

由于您的矩阵是对称的,因此 LU 分解实际上相当于 LDL* 分解,这就是您在实际情况中应该使用的分解。由于您的矩阵也是正定的,因此您应该能够在不进行旋转的情况下执行此分解。


编辑:请注意,对于此应用程序,您不需要解决完整的

C x = b
问题。

相反,给定

C = L D L*
和差向量
v = a-b
,求解
L* y = v
得到
y
(这是完整 LU 求解器工作量的一半)。

然后,可以在线性时间内计算

y^t D^-1 y =  v^t C^-1 v


11
投票

第一马哈拉诺比斯距离 (MD) 是相对于两个向量测量的不确定性的标准化距离。当

C=Indentity matrix
时,MD 减少到欧氏距离,因此乘积减少到向量范数。此外,对于所有非零向量,MD 始终是正定的或大于零。通过适当选择向量
a
b
进行公式化,
a*C^-1*b
可以小于零。希望您正在寻找的向量的差异是
x=(a-b)
,这使得计算
x^t*C^-1*x
,其中
x^t
是向量
x
的转置。另请注意
MD=sqrt(x^t*C^-1*x)
由于您的矩阵是对称且正定的,因此您可以利用 Cholesky 分解
(MatLab-chol)
,它使用一半的运算作为 LU,并且在数值上更稳定。
chol(C)=L
,其中
C=L*L^t
,其中
L
是下三角矩阵,
L^t
L
的转置,使其成为上三角矩阵。你的算法应该是这样的

(Matlab)

  x=a-b;
  L=chol(C); 
  z=L\x; 
  MD=z'*z; 
  MD=sqrt(MD);

1
投票
# Mahalanobis Distance Matrix

import numpy as np
from scipy.spatial import distance
from scipy.spatial.distance import mahalanobis
from scipy.spatial.distance import pdist
from scipy.spatial.distance import squareform

# Example
Data = np.array([[1,2],[3,2],[4,3]])
Cov = np.cov(np.transpose([[1,2],[3,2],[4,3]]))
invCov = np.linalg.inv(Cov)

Y = pdist(Data, 'mahalanobis', invCov)
MD = squareform(Y) 

0
投票

对于居中数据

import numpy as np
import pandas as pd

# Example
x = np.array([[1,2],[3,2],[4,3]])
df= df = pd.DataFrame(x,columns=['score',  'grade'])
print(df)

# first, to center data (or demean)
x_mu = x - np.mean(x)

Cov = np.cov(df.values.T)
invCov = np.linalg.inv(Cov)
left = np.dot(x_mu, invCov)
mahal = np.dot(left, x_mu.T)
m=  mahal.diagonal()
print(m)
© www.soinside.com 2019 - 2024. All rights reserved.