冗長性のあるワンホットエンコーディングによるパラメータ化の係数をフルランクの計画行列の係数に変換する

overparametrize.R

library(moltenNMF)
library(Matrix)

df = expand.grid(letters[1:2], LETTERS[1:2])
X = sparse_onehot(~., data = df)

A = matrix(
  c(0.5, 0.5, 0.5, 0.5, 0, 1, 0, 0, 0, 0, 0, 1),
  nrow = 4,
  ncol = 3
)
print(A)

Xs = sparse.model.matrix(~., data = df)
print(all(X %*% A == Xs))

B = matrix(
  c(1, -1, 0, 0, 1, 0, 1, 0, -1, 0, 0, 1),
  nrow = 3,
  ncol = 4,
)

print(B)
print(all(Xs %*% B == X))

df = expand.grid(letters[1:3], LETTERS[1:2])
Xs = sparse.model.matrix(~., data = df)
X = sparse_onehot(~., data = df)
B0 = matrix(
  c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0),
  nrow = 4,
  ncol = 5,
)

B1 = matrix(
  c(0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1),
  nrow = 4,
  ncol = 5,
)
B2 = matrix(
  c(0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0),
  nrow = 4,
  ncol = 5,
)


B0 + B1 - B2
clen = c(3, 3)

full_to_onehot <- function(clen) {
  n_full = sum(clen - 1) + 1
  n_onehot = sum(clen)
  lshift = 0
  rshift = 0
  A1 = matrix(0, n_full, n_onehot)
  for (k in seq_along(clen)) {
    A1[1, rshift + 1] = 1
    for (i in 2:clen[k]) {
      A1[lshift + i, rshift + i] = 1
    }
    lshift = lshift + clen[k] - 1
    rshift = rshift + clen[k]
  }

  A2 = matrix(0, n_full, n_onehot)
  lshift = 0
  rshift = 0
  for (k in seq_along(clen)) {
    for (i in 2:clen[k]) {
      A2[lshift + i, rshift + 1] = -1
    }
    lshift = lshift + clen[k] - 1
    rshift = rshift + clen[k]
  }
  A = A1 + A2
  A[1, 1] = 1
  return(A)
}

df = expand.grid(letters[1:4], LETTERS[1:2], factor(1:2))
Xs = sparse.model.matrix(~., data = df)
A = full_to_onehot(c(4, 2, 2))
Xs %*% A
###

set_attr_modelmat <- function(X) {
  if (!is.null(attr(X, "assign"))) {
    attr(X, "indices") <- c(0L, cumsum(rle(attr(X, "assign"))$lengths))
  }
  labs <- names(attr(X, "contrasts"))
  if (!is.null(labs)) {
    attr(X, "term.labels") <- labs
    if (!is.null(X@Dimnames[[2]])) {
      attr(X, "value.labels") <- sub(
        paste(labs, collapse = "|"),
        "",
        X@Dimnames[[2]]
      )
    }
  }
  return(X)
}

df = as.data.frame(Titanic)
f = Freq ~ Class + Sex + Age + Survived
X_s = sparse_onehot(f, data = df) #冗長なワンホット

L = 2
system.time({
  out_s <- mNMF_vb.default(df$Freq, X = X_s, L = L, iter = 1000)
})


X = sparse.model.matrix(f, data = df) #フルランク
X = set_attr_modelmat(X)

system.time({
  out <- mNMF_vb.default(df$Freq, X = X, L = L, iter = 1000)
})

plot(out_s$ELBO[-1], type = "l", lty = 2)
lines(out$ELBO[-1], type = "l")

V = (out$shape / out$rate)
V_s = (out_s$shape / out_s$rate)

fit = product_m.default(X, V)
fit_s = product_m.default(X_s, V_s)


abline(0, 1, lty = 2)


A = full_to_onehot(c(4, 2, 2, 2))
dim(A)
dim(V)

all(X %*% A == X_s)

V2 <- exp(A %*% log(V_s))
fit2 = product_m.default(X, V2)
plot(df$Freq, fit)
points(df$Freq, fit_s, col = "royalblue", pch = 2)
points(df$Freq, fit2, col = "orangered", pch = 3)

V = V[, order(colMeans(V), decreasing = TRUE)]
V2 = V2[, order(colMeans(V2), decreasing = TRUE)]
plot(V, V2)
abline(0, 1, lty = 2)

状況設定：冗長性のあるワンホットエンコーディングによるパラメータ化の係数をフルランクの計画行列の係数に変換したい．

ここで冗長性のあるワンホットエンコーディングは以下のようなことを指す．

カテゴリ変数が $K$ 水準あるとし，$K$ 個すべてのダミー変数を用いると，説明変数と回帰係数の一次結合は次のようになる．

$$ f = \sum_{k=1}^K \beta_k \mathbb{I}(x = k) \tag{1} $$

$\mathbb{I}(x = k)$ は指示関数（indicator function）とした．すなわち，$x = k$ なら 1，そうでなければ 0 の値を取る．

説明変数と回帰係数の一次結合 $f$ を 線形予測子 （linear predictor）と呼ぶ．

このとき，次が成り立つ．

$$ \sum_{k=1}^K \mathbb{I}(x=k)=1. $$

これは計画行列 $X$ の $(n,k)$ 成分を $x_{nk}=\mathbb{I}(x_n=k)$ とした場合，列基本変形ですべて0にできる列が存在することを意味する．

行列の列の数は $K$ だがランクは $K-1$ なので，このような計画行列を ランク落ち しているという．

ランク落ちしていない計画行列を フルランク であるという．

例えば，切片項 $\beta_0$ を取り，最初のカテゴリを基準にすると線形予測子は次のように書ける．

$$ f = \tilde{\beta}_0 + \sum_{k=2}^{K} \tilde{\beta}_k \mathbb{I}(x = k) \tag{2} $$

このパラメータの置き方に対応するフルランクの計画行列 $\tilde X$ は，$X$ からカテゴリ変数の2個め以降のカテゴリに相当する列を取ってくれば得られる．このことは $\tilde X = XA$ となる行列 $A$ があることを意味する．

いま，線形予測子 $f$ の値は(1)と(2)で完全に一致しているとする．

$$ \tilde {X} \tilde\beta = X \beta $$

である．

$$ XA \tilde\beta = X \beta $$

より，フルランク表現 $\tilde \beta$ から冗長な表現 $\beta$ への変換は

$$ \beta = A \tilde\beta $$

で得られる．

逆に， $X=\tilde XB$ となる行列 $B$ があれば，

$$ XB \tilde\beta = \tilde{X} \tilde{\beta} $$

より，フルランク表現 $\tilde \beta$ から冗長な表現 $\beta$ への変換は

$$ \tilde{\beta} = B \beta $$

で得られる．