如何计算非方阵的Cholesky分解,用numpy
计算马氏距离
def get_fitting_function(G):
print(G.shape) #(14L, 11L) --> 14 samples of dimension 11
g_mu = G.mean(axis=0)
#Cholesky decomposition uses half of the operations as LU
#and is numerically more stable.
L = np.linalg.cholesky(G)
def fitting_function(g):
x = g - g_mu
z = np.linalg.solve(L, x)
#Mahalanobis Distance
MD = z.T*z
return math.sqrt(MD)
return fitting_function
C:\Users\Matthias\CV\src\fitting_function.py in get_fitting_function(G)
22 #Cholesky decomposition uses half of the operations as LU
23 #and is numerically more stable.
---> 24 L = np.linalg.cholesky(G)
25
26 def fitting_function(g):
C:\Users\Matthias\AppData\Local\Enthought\Canopy\User\lib\site-packages\numpy\linalg\linalg.pyc in cholesky(a)
598 a, wrap = _makearray(a)
599 _assertRankAtLeast2(a)
--> 600 _assertNdSquareness(a)
601 t, result_t = _commonType(a)
602 signature = 'D->D' if isComplexType(t) else 'd->d'
C:\Users\Matthias\AppData\Local\Enthought\Canopy\User\lib\site-packages\numpy\linalg\linalg.pyc in _assertNdSquareness(*arrays)
210 for a in arrays:
211 if max(a.shape[-2:]) != min(a.shape[-2:]):
--> 212 raise LinAlgError('Last 2 dimensions of the array must be square')
213
214 def _assertFinite(*arrays):
LinAlgError: Last 2 dimensions of the array must be square
LinAlgError: Last 2 dimensions of the array must be square
基于Matlab的:协方差矩阵的Mahalanobis距离反演实现
编辑:chol(a)
= linalg.cholesky(a).T
cholesky分解矩阵( matlab中的chol(a)
返回上三角矩阵,linalg.cholesky(a)
返回下三角矩阵)(来源:链接)
Edit2:
G -= G.mean(axis=0)[None, :]
C = (np.dot(G, G.T) / float(G.shape[0]))
#Cholesky decomposition uses half of the operations as LU
#and is numerically more stable.
L = np.linalg.cholesky(C).T
所以如果D=x^t.S^-1.x=x^t.(L.L^t)^-1.x=x^t.L.L^t.x=z^t.z
发布于 2014-06-26 10:30:27
我不相信你能做到。Cholesky分解不仅需要一个方阵,而且需要一个Hermitian矩阵和一个正定矩阵来表示唯一性。它基本上是一个LU分解,但条件是L= U‘。事实上,该算法常被用作数值检验给定矩阵是否为正定的一种方法。见维基百科。
也就是说,从定义上说,协方差矩阵是对称正半定的,所以你应该能够在它上做cholesky。
编辑:当你计算的时候,你的矩阵C=np.dot(G, G.T)
应该是对称的,但是可能有什么不对。您可以尝试将其强制符号化为C = ( C + C.T) /2.0
,然后再试一次chol(C)
。
https://stackoverflow.com/questions/23685080
复制