如何计算非方阵的Cholesky分解，以便用“numpy”计算马氏距离

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:UsersMatthiasCVsrcfitting_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:UsersMatthiasAppDataLocalEnthoughtCanopyUserlibsite-packagesnumpylinalglinalg.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:UsersMatthiasAppDataLocalEnthoughtCanopyUserlibsite-packagesnumpylinalglinalg.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矩阵的choolesky分解(matlab中的chol(a)返回上三角矩阵，但linalg.cholesky(a)返回下三角矩阵)(来源:Link)

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 ^ -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)。