http://antkillerfarm.github.io/

因子分析的EM估计（续）

去掉和各参数无关的部分后，可得：

<script id="MathJax-Element-1" type="math/tex; mode=display">\begin{align} &\sum_{i=1}^mE\left[\log p(x^{(i)}\vert z^{(i)};\mu,\Lambda,\Psi)\right] \\&=\sum_{i=1}^mE\left[\frac{1}{(2\pi)^{n/2}\lvert\Psi\rvert^{1/2}}\exp\left(-\frac{1}{2}(x^{(i)}-\mu-\Lambda z^{(i)})^T\Psi^{-1}(x^{(i)}-\mu-\Lambda z^{(i)})\right)\right] \\&=\sum_{i=1}^mE\left[-\frac{1}{2}\log\lvert\Psi\rvert-\frac{n}{2}\log(2\pi)-\frac{1}{2}(x^{(i)}-\mu-\Lambda z^{(i)})^T\Psi^{-1}(x^{(i)}-\mu-\Lambda z^{(i)})\right] \end{align}</script>

去掉和 Λ <script id="MathJax-Element-2" type="math/tex">\Lambda</script>无关的部分，并求导可得：

<script id="MathJax-Element-3" type="math/tex; mode=display">\nabla_\Lambda\sum_{i=1}^m-E\left[\frac{1}{2}(x^{(i)}-\mu-\Lambda z^{(i)})^T\Psi^{-1}(x^{(i)}-\mu-\Lambda z^{(i)})\right]\tag{1}</script>

因为公式1中 E[⋅] <script id="MathJax-Element-4" type="math/tex">E[\cdot]</script>部分的结果实际上是个实数，因此该公式可变形为：

<script id="MathJax-Element-5" type="math/tex; mode=display">\nabla_\Lambda\sum_{i=1}^m-E\left[\operatorname{tr}\left(\frac{1}{2}(x^{(i)}-\mu-\Lambda z^{(i)})^T\Psi^{-1}(x^{(i)}-\mu-\Lambda z^{(i)})\right)\right]</script>

而：

<script id="MathJax-Element-6" type="math/tex; mode=display">\begin{align} &\frac{1}{2}(x^{(i)}-\mu-\Lambda z^{(i)})^T\Psi^{-1}(x^{(i)}-\mu-\Lambda z^{(i)}) \\&=\frac{1}{2}\left[((x^{(i)}-\mu)^T-(\Lambda z^{(i)})^T)\Psi^{-1}((x^{(i)}-\mu)-\Lambda z^{(i)})\right] \\&=\frac{1}{2}\left[(x^{(i)}-\mu)^T\Psi^{-1}(x^{(i)}-\mu)-(x^{(i)}-\mu)^T\Psi^{-1}\Lambda z^{(i)} \\-(\Lambda z^{(i)})^T\Psi^{-1}(x^{(i)}-\mu)+(\Lambda z^{(i)})^T\Psi^{-1}\Lambda z^{(i)}\right] \end{align}</script>

去掉和 Λ <script id="MathJax-Element-7" type="math/tex">\Lambda</script>无关的部分，可得：

<script id="MathJax-Element-8" type="math/tex; mode=display">\frac{1}{2}\left[(\Lambda z^{(i)})^T\Psi^{-1}\Lambda z^{(i)}-(x^{(i)}-\mu)^T\Psi^{-1}\Lambda z^{(i)}-(\Lambda z^{(i)})^T\Psi^{-1}(x^{(i)}-\mu)\right]</script>

所以：

<script id="MathJax-Element-9" type="math/tex; mode=display">\begin{align} &\nabla_\Lambda\sum_{i=1}^m-E\left[\operatorname{tr}\left(\frac{1}{2}\left[(\Lambda z^{(i)})^T\Psi^{-1}\Lambda z^{(i)}-(x^{(i)}-\mu)^T\Psi^{-1}\Lambda z^{(i)}-(\Lambda z^{(i)})^T\Psi^{-1}(x^{(i)}-\mu)\right]\right)\right] \\&\begin{split}=\nabla_\Lambda\sum_{i=1}^m-E\left[\frac{1}{2}\operatorname{tr}\left((\Lambda z^{(i)})^T\Psi^{-1}\Lambda z^{(i)}\right)-\frac{1}{2}\operatorname{tr}\left((x^{(i)}-\mu)^T\Psi^{-1}\Lambda z^{(i)}\right) \\-\frac{1}{2}\operatorname{tr}\left((\Lambda z^{(i)})^T\Psi^{-1}(x^{(i)}-\mu)\right)\right]\end{split} \\&=\sum_{i=1}^m\nabla_\Lambda E\left[-\frac{1}{2}\operatorname{tr}\left(\Lambda^T \Psi^{-1}\Lambda z^{(i)}(z^{(i)})^T\right)+\operatorname{tr}\left(\Lambda^T \Psi^{-1}(x^{(i)}-\mu)(z^{(i)})^T\right)\right]\tag{2} \end{align}</script>

因为：

<script id="MathJax-Element-10" type="math/tex; mode=display">\nabla_A\operatorname{tr}ABA^TC=CAB+C^TAB^T</script>

所以：

<script id="MathJax-Element-11" type="math/tex; mode=display">\begin{align} &\nabla_A\operatorname{tr}\left(\Lambda^T \Psi^{-1}\Lambda z^{(i)}(z^{(i)})^T\right)=z^{(i)}(z^{(i)})^T\Lambda^T \Psi^{-1}+(z^{(i)}(z^{(i)})^T)^T\Lambda^T(\Psi^{-1})^T \\&=z^{(i)}(z^{(i)})^T\Lambda^T \Psi^{-1}+((z^{(i)})^T)^T(z^{(i)})^T\Lambda^T\Psi^{-1}=2z^{(i)}(z^{(i)})^T\Lambda^T \Psi^{-1} \end{align}</script>

代入公式2，可得：

<script id="MathJax-Element-12" type="math/tex; mode=display">\sum_{i=1}^mE\left[-\Psi^{-1}\Lambda z^{(i)}(z^{(i)})^T+\Psi^{-1}(x^{(i)}-\mu)(z^{(i)})^T\right]</script>

由上式等于0，可得：

<script id="MathJax-Element-13" type="math/tex; mode=display">\sum_{i=1}^m\Lambda E_{z^{(i)}\sim Q_i}\left[z^{(i)}(z^{(i)})^T\right]=\sum_{i=1}^m(x^{(i)}-\mu)E_{z^{(i)}\sim Q_i}\left[(z^{(i)})^T\right]</script>

因此：

<script id="MathJax-Element-14" type="math/tex; mode=display">\Lambda=\left(\sum_{i=1}^m(x^{(i)}-\mu)E_{z^{(i)}\sim Q_i}\left[(z^{(i)})^T\right]\right)\left(\sum_{i=1}^m E_{z^{(i)}\sim Q_i}\left[z^{(i)}(z^{(i)})^T\right]\right)^{-1}\tag{3}</script>

因为：

<script id="MathJax-Element-15" type="math/tex; mode=display">\operatorname{Cov}(Y)=E[YY^T]-E[Y]E[Y]^T</script>

所以：

<script id="MathJax-Element-16" type="math/tex; mode=display">E[YY^T]=E[Y]E[Y]^T+\operatorname{Cov}(Y)</script>

因此根据之前的讨论可得：

<script id="MathJax-Element-17" type="math/tex; mode=display">E_{z^{(i)}\sim Q_i}\left[(z^{(i)})^T\right]=\mu_{z^{(i)}\vert x^{(i)}}^T</script>

<script id="MathJax-Element-18" type="math/tex; mode=display">E_{z^{(i)}\sim Q_i}\left[z^{(i)}(z^{(i)})^T\right]=\mu_{z^{(i)}\vert x^{(i)}}\mu_{z^{(i)}\vert x^{(i)}}^T+\Sigma_{z^{(i)}\vert x^{(i)}}</script>

将上式代入公式3，可得：

<script id="MathJax-Element-19" type="math/tex; mode=display">\Lambda=\left(\sum_{i=1}^m(x^{(i)}-\mu)\mu_{z^{(i)}\vert x^{(i)}}^T\right)\left(\sum_{i=1}^m \left(\mu_{z^{(i)}\vert x^{(i)}}\mu_{z^{(i)}\vert x^{(i)}}^T+\Sigma_{z^{(i)}\vert x^{(i)}}\right)\right)^{-1}</script>

这里需要注意的是，和之前的混合高斯模型相比，我们不仅要计算 Σz(i)|x(i) <script id="MathJax-Element-20" type="math/tex">\Sigma_{z^{(i)}\vert x^{(i)}}</script>，还要计算 E[z] <script id="MathJax-Element-21" type="math/tex">E[z]</script>和 E[zzT] <script id="MathJax-Element-22" type="math/tex">E[zz^T]</script>。

此外，我们还可得出：（推导过程略）

<script id="MathJax-Element-23" type="math/tex; mode=display">\mu=\frac{1}{m}\sum_{i=1}^mx^{(i)}</script>

<script id="MathJax-Element-24" type="math/tex; mode=display">\begin{split}\Phi=\frac{1}{m}\sum_{i=1}^m\left(x^{(i)}(x^{(i)})^T-x^{(i)}\mu_{z^{(i)}\vert x^{(i)}}^T\Lambda^T-\Lambda\mu_{z^{(i)}\vert x^{(i)}}(x^{(i)})^T \\+\Lambda\left(\mu_{z^{(i)}\vert x^{(i)}}\mu_{z^{(i)}\vert x^{(i)}}^T+\Sigma_{z^{(i)}\vert x^{(i)}}\right)\Lambda^T\right)\end{split}</script>

机器学习中的矩阵方法

在继续Andrew Ng的讲义之前，我们需要加强一些矩阵的相关知识。虽然Andrew Ng的讲义中已经包含了一个线性代数方面的简介文章，然而真的就只是简介而已，好多内容都没有。

这里推荐一本书《Matrix Methods in Data Mining and Pattern Recognition》。

作者：Lars Eld´en，执教于Linköping University数学系。

http://www.cnblogs.com/daniel-D/p/3204508.html

这是daniel-D写的中文笔记。

这一部分的内容属于数值计算领域，涉及的概念虽然不复杂，但提出一个高效算法，仍然不是件容易的事情。

还有另外一本书《Liner Algebra Done Right》，也值得推荐。这本书从定义矩阵算子，而不是通过行列式，来解释各种线性代数原理，提供了一种独特的视角。因为算子是有明确的几何或物理意义的，而行列式则不然。

作者：Sheldon Jay Axler，1949年生，美国数学家。普林斯顿大学本科，UCB博士，MIT博士后，San Francisco State University教授。美国的数学系基本就是本科和博士，很少有硕士。因为数学，尤其是理论数学，需要高度的抽象思维能力，半调子的硕士，既不好找工作，也不好搞科研。

三角矩阵的求逆问题

<script id="MathJax-Element-25" type="math/tex; mode=display">\begin{bmatrix} l_{11} & 0 & 0 \\ l_{21} & l_{22} & 0 \\ l_{31} & l_{32} & l_{33} \\ \end{bmatrix} \begin{bmatrix} u_{11} & u_{12} & u_{13} \\ 0 & u_{22} & u_{23} \\ 0 & 0 & u_{33} \\ \end{bmatrix} </script>

以3阶方阵为例，上面左边的矩阵被称为下三角矩阵（lower triangular matrix），而右边的矩阵被称为上三角矩阵（upper triangular matrix）。

对于矩阵求逆问题来说，下三角矩阵是一类比较简单的矩阵，求逆难度仅高于对角阵。

下三角矩阵的逆矩阵也是下三角矩阵，因此：

<script id="MathJax-Element-26" type="math/tex; mode=display">AA^{-1}=\begin{bmatrix} a_{11} & 0 & \dots & 0 \\ a_{21} & a_{22} & \dots & 0 \\ \dots & \dots & \dots & \dots \\ a_{n1} & a_{n2} & \dots & a_{nn} \\ \end{bmatrix} \begin{bmatrix} b_{11} & 0 & \dots & 0 \\ b_{21} & b_{22} & \dots & 0 \\ \dots & \dots & \dots & \dots \\ b_{n1} & b_{n2} & \dots & b_{nn} \\ \end{bmatrix} =\begin{bmatrix} 1 & 0 & \dots & 0 \\ 0 & 1 & \dots & 0 \\ \dots & \dots & \dots & \dots \\ 0 & 0 & \dots & 1 \\ \end{bmatrix} </script>

由矩阵乘法定义，可得：

<script id="MathJax-Element-27" type="math/tex; mode=display">c_{ij}=\sum_{k=j}^ia_{ik}b_{kj}</script>

由 cij=1,i=j <script id="MathJax-Element-28" type="math/tex">c_{ij}=1,i=j</script>，可得： bii=1aii <script id="MathJax-Element-29" type="math/tex">b_{ii}=\frac{1}{a_{ii}}</script>

由 cij=0,i≠j <script id="MathJax-Element-30" type="math/tex">c_{ij}=0,i\neq j</script>，可得：

<script id="MathJax-Element-31" type="math/tex; mode=display">c_{ij}=\sum_{k=j}^{i-1}a_{ik}b_{kj}+a_{ii}b_{ij}=0</script>

因此：

<script id="MathJax-Element-32" type="math/tex; mode=display">b_{ij}=-\frac{1}{a_{ii}}\sum_{k=j}^{i-1}a_{ik}b_{kj}=-b_{ii}\sum_{k=j}^{i-1}a_{ik}b_{kj}</script>

上三角矩阵求逆，可通过转置转换成下三角矩阵求逆。这里会用到以下性质:

<script id="MathJax-Element-33" type="math/tex; mode=display">(A^T)^{-1}=(A^{-1})^T</script>

LU分解

LU分解可将矩阵A分解为 A=LU <script id="MathJax-Element-34" type="math/tex">A=LU</script>，其中L是下三角矩阵，U是上三角矩阵。

LU分解的用途很多，其中之一是求逆：

<script id="MathJax-Element-35" type="math/tex; mode=display">A^{-1}=(LU)^{-1}=U^{-1}L^{-1}</script>

LU分解有若干种算法，常见的包括Doolittle、Cholesky、Crout算法。

注：Myrick Hascall Doolittlee，1830~1913。

Andr´e-Louis Cholesky，1875~1918，法国数学家、工程师、军官。死于一战战场。

Prescott Durand Crout，1907~1984，美国数学家，22岁获MIT博士。

这里只介绍一下Doolittle算法。

<script id="MathJax-Element-36" type="math/tex; mode=display">A=\begin{bmatrix} a_{11} & a_{12} & \dots & a_{1n} \\ a_{21} & a_{22} & \dots & a_{2n} \\ \dots & \dots & \dots & \dots \\ a_{n1} & a_{n2} & \dots & a_{nn} \\ \end{bmatrix}=LU= \begin{bmatrix} 1 & 0 & \dots & 0 \\ l_{21} & 1 & \dots & 0 \\ \dots & \dots & \dots & \dots \\ l_{n1} & l_{n2} & \dots & 1 \\ \end{bmatrix} \begin{bmatrix} u_{11} & u_{12} & \dots & u_{1n} \\ 0 & u_{22} & \dots & u_{2n} \\ \dots & \dots & \dots & \dots \\ 0 & 0 & \dots & u_{nn} \\ \end{bmatrix} </script>

由矩阵乘法定义，可知：

<script id="MathJax-Element-37" type="math/tex; mode=display">a_{1j}=u_{1j},j=1,2,\dots,n</script>

<script id="MathJax-Element-38" type="math/tex; mode=display">a_{ij}=\begin{cases} \sum_{t=1}^jl_{it}u_{tj}, & j

因此：

u1j=a1j,j=1,2,…,n <script id="MathJax-Element-39" type="math/tex">u_{1j}=a_{1j},j=1,2,\dots,n</script>（U的第1行）
lj1=aj1/u11,j=1,2,…,n <script id="MathJax-Element-40" type="math/tex">l_{j1}=a_{j1}/u_{11},j=1,2,\dots,n</script>（L的第1列）
For i=2,3,…,n <script id="MathJax-Element-41" type="math/tex">i=2,3,\dots,n</script> do
uii=aii−∑i−1t=1litutj <script id="MathJax-Element-42" type="math/tex">u_{ii}=a_{ii}-\sum_{t=1}^{i-1}l_{it}u_{tj}</script>
uij=aij−∑i−1t=1litutj <script id="MathJax-Element-43" type="math/tex">u_{ij}=a_{ij}-\sum_{t=1}^{i-1}l_{it}u_{tj}</script> for j=i+1,…,n <script id="MathJax-Element-44" type="math/tex">j=i+1,\dots,n</script>（U的第i行）
lji=aji−∑i−1t=1ljtutiuii <script id="MathJax-Element-45" type="math/tex">l_{ji}=\frac{a_{ji}-\sum_{t=1}^{i-1}l_{jt}u_{ti}}{u_{ii}}</script> for j=i+1,…,n <script id="MathJax-Element-46" type="math/tex">j=i+1,\dots,n</script>（L的第i列）
End
unn=ann−∑n−1t=1lntutn <script id="MathJax-Element-47" type="math/tex">u_{nn}=a_{nn}-\sum_{t=1}^{n-1}l_{nt}u_{tn}</script>

参见：

http://www3.nd.edu/~zxu2/acms40390F11/Alg-LU-Crout.pdf

QR分解

任意实数方阵A，都能被分解为 A=QR <script id="MathJax-Element-59" type="math/tex">A=QR</script>。这里的Q为正交单位阵，即 QTQ=I <script id="MathJax-Element-60" type="math/tex">Q^TQ=I</script>。R是一个上三角矩阵。这种分解被称为QR分解。

QR分解也有若干种算法，常见的包括Gram–Schmidt、Householder和Givens算法。

注：Jørgen Pedersen Gram，1850～1916，丹麦数学家，在矩阵、数论、泛函等领域皆有贡献。他居然是被自行车撞死的…

Erhard Schmidt，1876～1959，德国数学家，哥廷根大学博士，柏林大学教授。David Hilbert的学生。20世纪数学界的几位超级大神之一。1933年前的哥廷根大学数学系，秒杀其他所有学校。所谓“一流的学生去哥廷根，智商欠费才去藤校”。

Alston Scott Householder，1904～1993，美国数学家，芝加哥大学博士，田纳西大学教授。ACM主席。

James Wallace Givens, Jr.，1910～1993，美国数学家，普林斯顿大学博士，西北大学教授。参与UNIVAC I机器项目（1951年），这是最早的商用计算机。

这里只介绍Gram–Schmidt算法，这个算法虽然名为Gram–Schmidt，然而拉普拉斯和柯西早就已经用过了。

首先介绍一下向量的投影运算的符号表示。

如上图所示，根据余弦定理和向量点乘的定义可得：

a⋅b=|a||b|cosθ

<script id="MathJax-Element-50" type="math/tex; mode=display">a\cdot b=|a||b|\cos \theta</script>

因此，向量a在向量b上的投影向量 a1 <script id="MathJax-Element-51" type="math/tex">a_1</script>，可表示为：

a1=|a|cosθb^=|a|a⋅b|a||b|b|b|=a⋅b|b|2b=a⋅bb⋅bb=⟨a,b⟩⟨b,b⟩b

<script id="MathJax-Element-52" type="math/tex; mode=display">a_1=|a|\cos \theta\hat b=|a|\frac{a\cdot b}{|a||b|}\frac{b}{|b|}=\frac{a\cdot b}{|b|^2}b=\frac{a\cdot b}{b\cdot b}b=\frac{\langle a,b\rangle}{\langle b,b\rangle}b</script>

特别的，当b为单位向量时：

a1=⟨a,b⟩(4)

<script id="MathJax-Element-53" type="math/tex; mode=display">a_1=\langle a,b\rangle\tag{4}</script>

我们定义投影符号如下：

projea=⟨e,a⟩⟨e,e⟩e

<script id="MathJax-Element-54" type="math/tex; mode=display">\mathrm{proj}_{\mathbf{e}}\mathbf{a} = \frac{\left\langle\mathbf{e},\mathbf{a}\right\rangle}{\left\langle\mathbf{e},\mathbf{e}\right\rangle}\mathbf{e}</script>