如何用`numpy`计算非方阵的Cholesky分解以求Mahalanobis距离？

Question

如何计算一个非方阵的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的实现，具体可以参考这个链接：马哈拉诺比斯距离与协方差矩阵的反演

补充说明：在chol(a)中，= linalg.cholesky(a).T，这是对一个矩阵进行Cholesky分解的方式（在Matlab中，chol(a)返回的是一个上三角矩阵，而linalg.cholesky(a)返回的是一个下三角矩阵）（来源：链接）

补充说明2：

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

Answer 1

我觉得你可能做不到。Cholesky分解不仅需要一个方阵，还需要它是厄米矩阵，并且要是正定矩阵才能保证结果的唯一性。简单来说，这个过程类似于LU分解，但有个条件是L等于U的转置。实际上，这个算法常常用来数值上检查一个矩阵是否是正定的。你可以查看一下维基百科。

不过，协方差矩阵根据定义是对称的正半定矩阵，所以你应该可以对它进行Cholesky分解。

补充一下：当你计算的时候，你的矩阵C=np.dot(G, G.T)应该是对称的，但可能有什么地方出错了。你可以尝试强行让它对称，方法是C = ( C + C.T) /2.0，然后再试试chol(C)。