[Julia] 行列分解を書いてみる(関係データ学習 その2)


最近久しぶりにやっていました.まだ問題は解決してない気もしますが,途中で一度まとめておきます.


4章の行列分解アルゴリズムのうち,一階の情報のみを使うアルゴリズム4.1を説明します.

(I, J)行列 \(X\) が与えられたときに、2つの行列; (I, R)行列\(U\) と (J, R)行列\(V\)を用意して,\(X\approx UV^T\)の形で近似します.RはIやJよりも小さい値に設定することで,情報をある意味で圧縮したことに相当します.

  • 行列のノルムにフロベニウスノルムを採用した上で,この問題は次のように表すことができます.

\[ \min_{U, V} ||X - UV^T||^2 _\mathrm{Fro} \]

  • 正則化項を導入する場合は次のように表すことができます.

\[ \min_{U, V} ||X - UV^T||^2 _\mathrm{Fro} + \lambda \left( ||U||^2 _\mathrm{Fro} + ||V||^2 _\mathrm{Fro}\right) \]

  • 具体的には行列U/Vを適当に初期化した後に,次のように計算します(アルゴリズム4.1).
1. U, Vを乱数で初期化
2. repeat
3.  g <- f(U, V)
4.  U <- U - eta * nablaU f(U, V)  [定義] -XV + U(V' V + λI)
5.  V <- V - eta * nablaV f(U, V)  [定義] -X'U + V(U'U + λI)
6. until | g - f(U, V) | / f(U, V) < ε

  • 実際に実装してみます.Juliaでは一部のギリシャ文字がそのまま使えるので,使っています.
  • 入力データはJLDファイルに書き込まれているとしています.

Juliaでの実装

using JLD

X = jldopen("data/sample_org.jld", "r") do file
  read(file, "Xorg")
end

N, M = size(X)
η = 0.0000001 # learning parameter
λ = 0.1 # regularization constant
K = 20
# K次元埋め込み
U = rand(N, K) * 50
V = rand(M, K) * 50

function Fro(X)
  sum = 0.0
  N, M = size(X)
  for i=1:N, j=1:M
    sum += (X[i, j] ^ 2)
  end
  return sqrt(sum)
end

# 目的関数
function fl2(X, U, V, λ)
  N, M = size(X)
  ret = 0.0
  counter = 0
  for i=1:N, j=1:M
    ret = (X[i, j] - U[i,:] ⋅ transpose(V[j, :])) ^ 2
  end
  ret *= 0.5

  # regularization term
  value = 0.0
  for k=1:K
    for i = 1:N
      value += U[i, k]^2
    end
    for j = 1:M
      value += V[j, k]^2
    end
  end
  ret += λ * 0.5 * value
  return ret
end

# 微分
function ∇u(X, U, V, λ, K)
  return -X * V + U * (transpose(V) * V + λ .* eye(K))
end

function ∇v(X, U, V, λ, K)
  return -transpose(X) * U + V * (transpose(U) * U + λ .* eye(K))
end

# update until convergence
n_iter = 1
while true
  g = fl2(X, U, V, λ)
  # println("Fro $(Fro(U))")
  U -= η .* ∇u(X, U, V, λ, K)
  V -= η .* ∇v(X, U, V, λ, K)
  gg = fl2(X, U, V, λ)
  println("$n_iter $g -> $gg")
  if abs(g - gg) / gg < 1e-3
    break
  end
  n_iter += 1
end
  • コードをそのまま書きました.
  • 学習率を0.0000001などの小さい値にしないと速攻で発散して死にます
  • 正則化パラメータを0.1, 0.01, 0.0で画像出力してみました.

0.1の場合

0.01の場合

0の場合

  • 結論: よく分からない
    • 本などを読んでも簡単に動くようなイメージがあるが,発散してしまいイマイチ実装がうまくいかない.
    • よくある二乗誤差でMSEを使う場合だと,結構うまくいく.ただし正則化をかけると死ぬことがあった.
  • 引き続き調査していきます.
  • 実装の比較や,よくあるMovieLenz,テンソル版などを2月〜3月でまとめる予定.