正如在高斯混合模型中推导的结果一样，我们在计算对数似然值log-likelihood和计算responsibility即gamma(z_nk)时，需要计算高斯分布的概率密度。并且正如推导过程一样，我们可以计算精度矩阵的行列式的倒数，而非直接计算协方差矩阵的逆阵。
我们先来计算密度函数的指数部分，这部分仅仅计算一个二次型在一个向量上的函数值（一个双层循环），再对其取指数函数得到。当然，考虑到精度矩阵的对称性，我们只需计算对角线部分以及对角线上三角部分即可。
import math
def exponent_of_gaussian(x, mean, precision):
x = x - mean
dim = x.shape[0]
quad = 0.0
# 上三角部分（不包括对角线）
for i in range(0, dim-1): # 从第一行到倒数第二行
for j in range(x.shape[0]): # 从第二列到最后一列
quad += precision[i][j] * x[i] * x[j]
# 乘以2得到二次型中除主对角线以外的求和项
quad *= 2
# 对角线部分
for i in range(dim):
quad += precision[i][i] * x[i] * x[i]
return math.exp(-0.5 * quad)
- 我们再来计算密度函数的系数部分，语音识别的论文称为G-Const，意思大概就是Gauss分布的Constant部分？
import numpy as np
def constant_of_gaussian(mean, precision):
dim = mean.shape[0]
g_const = math.pow(2*math.pi, dim)
g_const /= np.linalg.det(precision)
return math.pow(g_const, -0.5)
- 当然，这只是我随便调了几个python函数计算出来的，实际计算要考虑计算复杂度，一般会用C语言实现。还要避免冗余计算和考虑数值稳定性。其中最主要的部分就是计算精度矩阵的行列式，可以使用矩阵分解算法LU、Cholesky来计算。
- LU大概简述一下如下：对任意一个矩阵A，可由高斯消元法进行化简，即存在一组基础矩阵P1,P2,...,PN使得PN...P2P1A的结果是reduced row echelon form矩阵R，若一个方阵A是非退化的，则R是一个上三角阵U且对角线不为0；如果这些基础矩阵中不包括交换两行的矩阵，则L_inv=PN...P2*P1是一个下三角矩阵，而U是上三角；又基础矩阵是可逆阵，所以L_inv可逆，两边乘以L_inv的逆矩阵L，就得到A=LU，易证一个可逆的下三角矩阵的逆矩阵也是下三角阵。
- 但是一般来说P1,P2,...,PN会包括交换两行的矩阵，而交换两行的矩阵不是下三角阵，因而对一般的矩阵（甚至满秩）来说A的LU分解不成立。因而问题就出在交换行这种基础矩阵，而我们不需要交换矩阵A任意两行时，就意味着在使用高斯消元法构造reduced row echelon form时，拿到的pivots均不为0。而这就意味着对任意阶主子矩阵A[0:k-1, 0:k-1]， k<n，它的行列式全不为0。
- 进一步，加强A的条件：当A是实对称矩阵正定矩阵时，A的所有pivots不为0，因此可以进行LU分解；不仅如此，还有U=L_tran（L的转置），即A=L*L_tran；这被称为Cholesky分解，理论上有着更强的条件，会有更快速的分解算法。
- 这里假设我们的协方差矩阵是实对称正定矩阵，则精度矩阵P也是实对称正定阵，P分解成LU后，下三角矩阵L的对角线全是1，根据行列式的性质，P的行列式就等于U的主对角线上的元素乘积。根据这位网友的分解图示可以容易写出示例程序：
def LU_decomposition(A):
if A.shape[0] != A.shape[1]:
raise DimensionError("must be a square matrix")
dim = A.shape[0]
L = np.eye(dim, dim) # L的主对角线全是1
U = np.zeros(dim, dim)
U[0, :] = A[0, :] # U的第一行就是A的第一行
L[1:, 0] = A[1:, 0] / U[0, 0] # L的第一列可用A的第一列除以A[0,0]得到
for i in range(1, dim): #从第二行开始，以行为主循环
# 观察图中的分解可知，虽然可以将行和列处理称统一的形式，但为了尽量减少除法运算，还是对每个i，将行和列分开处理
for j in range(i, dim): # 处理行，注意这里使用np.dot库函数也只是表达含义，而非效率，工程中肯定要用blas专业的函数
U[i, j] = A[i, j] - np.dot(L[i, 0:i], U[0:i, j])
for j in range(i+1, dim): # 处理列，需要额外一次除法运算
L[j, i] = (A[j, i] - np.dot(L[j, 0:i-1], U[0:i-1, i])) / U[i, i]
return L, U
def determinant(A):
_, U = LU_decomposition(A)
return np.product(np.diagonal(U))
- 回到LU算法，实际应用中不可能要求A有这么强的条件，这样会限制使用范围。我们在使用Gauss消元法过程中，一定可以用一系列基础矩阵左成A，将A化为上三角阵，只不过需要额外的交换行基础矩阵而已，假设在对A进行Gauss消元法过程中，涉及的所有交换行基础矩阵的乘积为P，事实上可以先对A左乘P进行换行，然后就可以进行上面基本的LU分解了，即：PA = LU，此时A为非奇异当且仅当U非奇异，且det(A) = s det(U)，s 是由permutation矩阵P确定的符号。
