MM Algorithm

Introduction
- MM 알고리즘을 왜 사용하는가?
Majorization과 Minorization
Example
- Least Square without Matrix Inversion
Median Regression
Reference

이 글은 고려대학교 통계학과 신승준 교수님의 2018년 1학기 통계 계산 방법론(대학원) 강의를 바탕으로 작성되었습니다.

Introduction

MM 알고리즘은 Majorization-Minimization, Minorization-Maximization의 약자이다.
MM 알고리즘은 high dimensional setting에 매우 적합하다.
Newton-Rhapson 같은 고전적인 optimization 방법들은 매우 큰 matrix를 다뤄야 하거나 matrix inverse 연산을 해야 하기 때문에 computation cost가 매우 크다. 이러한 점을 해결하기 위해 나온 방법 중 하나가 MM 알고리즘이다.

Majorization

Minorizaiton

MM 알고리즘을 왜 사용하는가?

1. It can generate an algorithm that avoids matrix inversion. + parameter를 분리해서 매우 큰 matrix에서 inverse 연산을 하는 상황을 피함

2. It can separate the parameters of a problem. + 한 번에 parameter 한 개씩 업데이트 + high dimensional problem을 일련의 one(또는 low) dimentional problem으로 바꿈

3. It can linearize an optimization problem. + derivative가 linear인 quadratic surrogate function을 이용

4. It can deal gracefully with equality and inequality constraints. + surrogate function의 optimization을 정의할 때 Jensen’s Inequality 같은 개념이 들어감

5. It can restore symmetry + 대부분의 closed form solution들은 symmetry에 의존함 + symmetry로 만들 수 있으면 surrogate function을 최적화하기 위한 closed form solution을 얻을 수 있음 + objective function이 symmetry가 아니어도 surrogate function을 symmetry로 만들면 closed form solution을 얻을 수 있음

6. It can turn a non-differentiable problem into a smooth problem. + objective function이 non-differentiable이어도 surrogate function을 differentiable(smooth)로 잡으면 미분의 여러가지 도구를 사용가능

Majorization과 Minorization

Majorization

\(x^{(t)}\)가 t번째 iteration에서 x값을 의미할 때, 다음 조건을 만족하면 함수 \(g(x|x^{(t)})\)가 함수 \(f(x)\)를 majorize한다고 말한다.

(i) \(g(x^{(t)}|x^{(t)}) = f(x^{(t)})\)

\(x=x^{(t)}\)일 때 g는 f에 닿는다.

(ii) \(g(x|x^{(t)}) \ge f(x^{(t)}) \quad \forall x\)

모든 x에 대해 g가 f보다 위에 있다.

[Figure1] Majorization

Figure1에서 g는 surrogate function이고 f는 objective function

Minorization

Minorization은 Majorization의 반대라고 생각하면 된다.

\(x^{(t)}\)가 t번째 iteration에서 x값을 의미할 때, 다음 조건을 만족하면 함수 \(g(x|x^{(t)})\)가 함수 \(f(x)\)를 minorize한다고 말한다.

(i) \(g(x^{(t)}|x^{(t)}) = f(x^{(t)})\) + \(x=x^{(t)}\)일 때 g는 f에 닿는다.

(ii) \(g(x|x^{(t)}) \le f(x^{(t)}) \quad \forall x\) + 모든 x에 대해 g가 f보다 아래에 있다.

Majorization-Minimization

\(g(x|x^{(t)})\)가 \(f(x)\)를 majorize하는 것은 \(-g(x|x^{(t)})\)가 \(-f(x)\)를 minorize하는 것으로 생각할 수 있다.

Majorization-Minimization은 \(g(x|x^{(t)})\)가 \(f(x)\)를 majorize할 때 다음 식을 푸는 것이다.

\[x^{(t+1)}=\underset{x}{argmin}\;g(x|x^{(t)})\]

실제 objective function인 \(f(x)\)대신 이를 Majorize하는 surrogate function \(g(x|x^{(x)})\)를 최소화한다.
Minorization-Maximization은 반대로 생각해서 argmax로 계산하면 된다.

Descent Property

\(x^{(t+1)}\)이 \(g(x|x^{(t)})\)를 최소화하는 지점에서의 x값이므로 다음의 관계가 성립한다는 것을 쉽게 알 수 있다.

\[f(x^{(t+1)}) \le f(x^{(t)})\]

proof

\(f(x^{(t+1)}) = g(x^{(t+1)}|x^{(t)})\)

\(\quad \quad \quad \quad \le g(x^{(t)}|x^{(t)}) =f(x^{(t)})\)

Example

Least Square without Matrix Inversion

convex fucntion \(f(t)\)는 다음과 같이 정의된다.

\[f(\alpha^Tt) \le \underset{j=1}{\overset{p}{\Sigma}}\alpha_jf(t_j)\]

이를 LS문제의 notation에 맞게 바꾸면

\[\underset{i=1}{\overset{n}{\Sigma}}f(x_i^T\beta)=\underset{i=1}{\overset{n}{\Sigma}}(y_i-x_i^T\beta)^2 \le \underset{i=1}{\overset{n}{\Sigma}} \underset{j=1}{\overset{p}{\Sigma}} \frac{x_j\beta_j^{(t)}}{x^T\beta^{(t)}}f(\frac{x^T\beta^{(t)}}{\beta_j^{(t)}}\beta_j)\]

이런 형태의 부등식을 만들 수 있다. 여기서 \(x\)와 \(\beta\)는 모두 양수여야 하는데 이 조건을 풀어주면

\[f(x^T\beta) \le \underset{j=1}{\overset{p}{\Sigma}}\alpha_j f(\frac{x_j}{\alpha_j}(\beta_j - \beta_j^{(t)})+ x^T\beta^{(t)})=g(\beta|\beta^{(t)})\]

where \(\alpha_j = \frac{|x_j|^p}{\underset{j=1}{\overset{p}{\Sigma}}|x_j|^p}\)

형태로 만들 수 있고,

이를 다시 쓰면

\[\underset{i=1}{\overset{n}{\Sigma}}f(x_i^T\beta)=\underset{i=1}{\overset{n}{\Sigma}}(y_i-x_i^T\beta)^2 \le \underset{i=1}{\overset{n}{\Sigma}}\underset{j=1}{\overset{p}{\Sigma}}\alpha_{ij} (y_i - \frac{x_{ij}}{\alpha_{ij}}(\beta_j - \beta_j^{(t)}) - x_i^T\beta^{(t)})^2\]

이 된다.

이제 \(\underset{i=1}{\overset{n}{\Sigma}}\underset{j=1}{\overset{p}{\Sigma}}\alpha_{ij} (y_i - \frac{x_{ij}}{\alpha_{ij}}(\beta_j - \beta_j^{(t)}) - x_i^T\beta^{(t)})^2\)를 minimize하면 된다.

위의 식을 미분해서 0이라고 놓고 풀면

\[\beta_j^{(t+1)}=\beta_j^{(t)} + \frac{\underset{i=1}{\overset{n}{\Sigma}}x_{ij}(y_i-x_i^T\beta^{(t)})}{\underset{i=1}{\overset{n}{\Sigma}}x_{ij}^2/\alpha_{ij}}\]

을 얻을 수 있다.

이제 \(\beta^{(t+1)}\)과 \(\beta^{(t)}\)의 차이가 매우 작아질 때까지 여러 번 iteration을 거치면 \(\hat{\beta}\)를 얻을 수 있다.

set.seed(1)
n <- 100
p <- 5

beta <- rep(1:p) # without intercept

x <- matrix(rnorm(n*p), n, p)
e <- rnorm(n, 0, 0.2)

y <- x %*% beta + e

mm.ls <- function(x, y, max.iter = 100, eps = 1.0e-8)
{
  beta <- rep(0, p)
  iter <- 1
  for (iter in 1:max.iter)
  {
    alpha <- abs(x)^1 / apply(abs(x)^1, 1, sum) # p=1???
    r <- c(y - x %*% beta)
    temp1 <- apply(x * r, 2, sum)
    temp2 <- apply(x^2/alpha, 2, sum)
    beta.new <- beta + temp1/temp2
    if (max(abs(beta.new - beta)) < eps) break
    beta <- beta.new
  }
  return(beta.new)
}

\(x : n \times p \;\; matrix\)
\(beta : p \times1 \;\; vector\)

# existing function
beta.ls <- coef(lm(y ~ -1 + x))
# mm algorithm
beta.ls.mm <- mm.ls(x, y) 
print(round(cbind(beta.ls, beta.ls.mm), 5))

##    beta.ls beta.ls.mm
## x1 0.97674    0.97674
## x2 1.99655    1.99655
## x3 2.99468    2.99468
## x4 3.98047    3.98047
## x5 4.98872    4.98872

alpha 계산할 때 p=1, 2는 같은 값이 나오는데 3부터 값이 많이 달라짐

Median Regression

Median regression의 objective function은 아래와 같다.

\[h(\beta)=\underset{i=1}{\overset{n}{\Sigma}}|y_i - x_i^T\beta|=\underset{i=1}{\overset{n}{\Sigma}}|r_i(\beta)|=\underset{i=1}{\overset{n}{\Sigma}}\sqrt{[r_i(\beta)]^2}\]

r_i() : i번째 residual

Majorization에 사용할 inequality는 아래와 같은데, 어디서 가져왔는지는 모르겠지만 아마 Approximation을 한 것이든지 아니면 유명한 inequality 중 하나일 것이다.

\[\sqrt{u} \le \sqrt{u^{(t)}} + \frac{u-u^{(t)}}{2\sqrt{u^{(t)}}}\]

위의 부등식을 이용한 Majorization은 아래와 같다.

\[h(\beta) \le \underset{i=1}{\overset{n}{\Sigma}}[\sqrt{[r_i(\beta^{(t)})]^2} + \frac{[r_i(\beta)]^2-[r_i(\beta^{(t)})]^2}{2\sqrt{[r_i(\beta^{(t)})]^2}}]\]

결과적으로 \(\beta\)는 다음 식을 iteratively 풀면 얻을 수 있다.

\[\beta^{(t+1)}=\underset{\beta}{argmin} \underset{i=1}{\overset{n}{\Sigma}} \frac{(y_i -x_i^T\beta)^2}{2\sqrt{[r_i(\beta^{(t)})]^2}}\]

\(2\sqrt{[r_i(\beta^{(t)})]^2}\)가 상수이므로 WLS 문제를 푸는 것과 같다.

mm.md <- function(x, y, max.iter = 100, eps = 1.0e-8)
{
  beta <- rep(0, p)
  iter <- 1
  for (iter in 1:max.iter)
  {
    r <- c(y - x %*% beta)
    w <- 1/abs(r + 1.0e-12) # avoid overflow => 0이면 inf가 나오는 경우가 있어서

    # QR decomposition 이용해서 WLS 푸는 것
    tilde.x <- x * sqrt(w)
    tilde.y <- y * sqrt(w)
    qr.obj <- qr(tilde.x)
    beta.new <- backsolve(qr.obj$qr, qr.qty(qr.obj, tilde.y))
    if (max(abs(beta.new - beta)) < eps) break
    beta <- beta.new
  }
  return(c(beta.new))
}

library(quantreg)

## Loading required package: SparseM

## 
## Attaching package: 'SparseM'

## The following object is masked from 'package:base':
## 
##     backsolve

# existing function
beta.md <- coef(rq(y ~ -1 + x))
# mm algorithm
beta.md.mm <- mm.md(x, y)

print(round(cbind(beta.md, beta.md.mm), 5))

##    beta.md beta.md.mm
## x1 0.98665    0.98699
## x2 2.01027    2.00402
## x3 2.95653    2.95774
## x4 3.96316    3.96799
## x5 5.01686    5.00995

rq()함수는 디폴트가 median regression이다.
값이 조금 차이가 나는 것은 numerical stability 때문이다.

Reference

Hunter, Lange(2003). A tutorial on MM algorithm.

Kenneth Lange의 강의 슬라이드

Kenneth Lange의 Youtube 강의