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

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), 4,3)

print(A)

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

###

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)

plot(df$Freq, fit)
points(df$Freq, fit_s, col="royalblue", pch=2)
abline(0,1,lty=2)

fullranktrans <- function(X){
  ind = attr(X, "indices")
  termlen = length(ind[-1])
  A = matrix(0,ind[length(ind)], ind[length(ind)]-termlen+1)
  for(k in 1:termlen){
    start = ind[k] + 2
    for(i in start:ind[-1][k]){
      A[i,i-k+1] <- 1
    }
  }
  A[,1] <- 1/termlen
  return(A)
}

A = fullranktrans(X_s)
print(all(X==X_s%*%A))
#[1] TRUE

V2 =exp(A%*%log(V))
image(V_s)
image(V2)
cor(V_s,V2)

fit2 = product_m.default(X_s, V2)
plot(df$Freq, fit2)
points(df$Freq, fit_s, col="royalblue", pch=2)
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 $$

で得られる．