前面介绍过EM算法，本文中介绍的MM算法是EM算法的一种更一般形式。当需要优化的问题很复杂，不能得到封闭形式的解时，一个好的数值解通常是需要的。MM算法将原优化问题转化为一系列的简单优化问题，使得问题的求解变得简单。MM的名称起着两个作用：对于最小化问题，分别代表Majorize和Minimize；对于最大化问题，分别代表Minorize和Maximize。后面介绍Majorize和Minorize的概念。
MM算法用一系列简单的优化问题代替一个困难的优化问题。简单体现在（1）避免大矩阵求逆，（2）线性化优化问题，（3）分离优化问题的参数，（4）更优雅的处理等式和不等式约束，（5）将一个不可微问题转化为一个光滑问题。迭代计算则是我们需要付出的代价。
MM算法思想：令$\theta^{(m)}$是参数$\theta$的一个固定值，$g(\theta|\theta^{(m)})$表示$\theta$的实值函数，其形式依赖于$\theta^{(m)}$。定义函数$g(\theta|\theta^{(m)})$在点$\theta^{(m)}$处Majorize函数$f(\theta)$，如果成立
\[g(\theta|\theta^{(m)})\ge{f(\theta)},\ \forall\ \theta,\]
\[g(\theta^{(m)}|\theta^{(m)})=f(\theta^{(m)}).(*)\]
即曲面$\theta\mapsto{g(\theta|\theta^{(m)})}$在曲面$f(\theta)$上方且相切于点$\theta=\theta^{(m)}$。定义函数$g(\theta|\theta^{(m)})$在点$\theta^{(m)}$处Minorize函数$f(\theta)$如果$-g(\theta|\theta^{(m)})$在点$\theta^{(m)}$处Majorize函数$-f(\theta)$.
通常$\theta^{(m)}$表示当前迭代开始时原优化问题的解。对于Majorize-Minimize MM算法，我们最小化majorizing函数$g(\theta|\theta^{(m)})$来代替最小化$f(\theta)$。令$\theta^{(m+1)}$为$g(\theta|\theta^{(m)})$的极小子，则MM算法使得$f(\theta)$沿着下山方向前进，即
\begin{align*}
f(\theta^{(m+1)})&=g(\theta^{(m+1)}|\theta^{(m)})+f(\theta^{(m+1)})-g(\theta^{(m+1)}|\theta^{(m)})\\
&\le{g(\theta^{(m)}|\theta^{(m)})+f(\theta^{(m)})-g(\theta^{(m)}|\theta^{(m)})}\\
&=f(\theta^{(m)})(**)
\end{align*}
这可由$g(\theta^{(m+1)}|\theta^{(m)})\le{g(\theta^{(m)}|\theta^{(m)})}$和(*)式得到。
下降性质(**)保证了MM算法的数值稳定性。对于最大化问题，可以类似处理，即最大化minorizing函数$g(\theta|\theta^{(m)})$得到下一次迭代的参数估计值$\theta^{(m+1)}$。
下面应用MM算法求解Logistic回归的参数估计。Logistic回归是一种响应变量为二元数据的广义线性模型。令$Y$为$n\times{1}$的响应变量，$X$为$n\times{p}$的自变量，Logistic模型假设$Y_i=1$的概率为
\[\pi_i(\theta)=exp\{x_i^T\theta\}/(1+exp\{x_i^T\theta\})。\]
则参数$\theta$的极大似然为
\begin{align*}
f(\theta)&=\prod\limits_{i=1}^n(\frac{e^{x_i^T\theta}}{1+e^{x_i^T\theta}})^{y_i}(\frac{1}{1+e^{x_i^T\theta}})^{1-y_i}\\
&=\prod\limits_{i=1}^n\frac{e^{x_i^T\theta{y_i}}}{1+e^{x_i^T\theta}}
\end{align*}
其对数极大似然为
\[l(\theta)=\sum\limits_{i=1}^n[Y_ix_i^T\theta-\log(1+e^{x_i^T\theta})]\]
我们要最大化$l(\theta)$。通常的Newton-Raphson算法迭代公式：
\[\theta_{m+1}=\theta_m-(\nabla^2l(\theta_m))^{-1}\nabla{l(\theta_m)}\]
对于每次新的$\theta_{m}$，我们需要重新计算矩阵$\nabla^2l(\theta_m)$的逆，而矩阵求逆是计算量大的操作。注意
\[\nabla{l(\theta)}=\sum\limits_{i=1}^n[Y_i-\pi_i(\theta)]x_i\]
\[\nabla^2l(\theta)=\sum\limits_{i=1}^n-\pi_i(\theta)[1-\pi_i(\theta)]x_ix_i^T\]
因为$\pi_i(\theta)[1-\pi_i(\theta)]$的上界是$1/4$，定义负定矩阵$B=-\frac{1}{4}X^TX$，则$\nabla^2l(\theta_m)-B$为一个非负定矩阵，定义
\[g(\theta|\theta^{(m)})=l(\theta^{(m)})+\nabla{l(\theta^{(m)})}^T(\theta-\theta^{(m)})+\frac{1}{2}(\theta-\theta^{(m)})^TB(\theta-\theta^{(m)})\]
则$g(\theta|\theta^{(m)})$是$l(\theta)$的minorizing函数，极大化$g(\theta|\theta^{(m)})$得到
\begin{align*}
\theta^{(m+1)}&=\theta^{(m)}-B^{-1}\nabla{l(\theta^{(m)})}\\
&=\theta^{(m)}-4(X^TX)^{-1}X^T[\pi(\theta^{(m)})-Y]
\end{align*}
其中$\pi(\theta_{(m)})=(\pi_1(\theta^{(m)}),\cdots,\pi_n(\theta^{(m)}))^T$.
应用MM算法，得到了一个更好的迭代公式，这个公式中$X^TX$的逆只需要计算一次即可，当$p$很大时，相对于Newton-Raphson算法优势明显。