Julia言語 Randomized eigenvalue decomposition を実装した．

n×nの実対称行列Aの固有値を大きい順にr<n個得たいです．lapackを使うと

val, vec = eigen(Symmetric(A), n-r+1:n)

でOKです．この時，eigenの第一引数はSymmetric型じゃないといけないので注意．（ところで，ここでeigenではどんなアルゴリズムが動いてるんですかね？全部の固有値を求める方法は調べると出てきますが，大きい順に，っていう時は何が裏で動いているんでしょう．．．？知っている人がいたら教えてほしいです．）

またArpackを使って，

using Arpack
val, vec=eigs(A, nev=r, which=:LR);

でもOKでした．ArpackはKrylov部分空間法の固有値解法ライブラリで，疎行列でも動くライブラリみたいです．このブログではArpackのことを何度か「疎行列向けのライブラリ」と表現しましたがそれはミスリーディングっぽいです．ごめんなさい．ここで，eigs(A)のAは実対称行列じゃなくても動きます．

自分のいま取り組んでいる研究で，Aの固有値全部はいらないから，大きい順にr個だけを爆速で取り出したい，という要請があります．こないだちょっと試してみたところ，eigen(Symmetric(A), n-r+1:n)とeigs(A, nev=r, which=:LR)は固有値を全部求めるeigen(A)と比べても，そんなに速くないみたいでした．そこで，噂に聞く結構速いアルゴリズムであるらしい Randomized eigenvalue decomposition(RED) を調べてみました．REDの特異値分解版のRSVDに関する実装や説明は沢山あるみたいですけど，REDに関しては，少なくともJulia実装は見当たりませんでした．一応，RandomizedLinAlg.jlというJuliaのライブラリがありますが，試してみても現在は少なくとも密行列についてはreigenは動かないです．（疎行列については分からない）

最適化済みLAPACKはSIMD、スレッド並列化が両方効いて相当速いが，ARPACKは並列化がうまく効かないから遅いらしいです．一方で，REDを使うと、ARPACKと計算量は同じですが、r≠1ではBLAS3を使える分、SIMDが効いて早くなることが期待されるらしいです．これは強い人に教わりました．

そこで次の2つの論文を参考に，REDを自分で実装してみました．固有値分解の周辺のトピックが全然分からないので論文読んだり調べたりしながら実装するのは，大変でした．１０数行のコードを書くのに４日以上かかりました．．．

epubs.siam.org

arxiv.org

ちなみに，１つ目のN.Halko論文に関連して，この方のD論はこちらです．

www.proquest.com

REDのJulia実装

見よう見まねで書いたRandomized eigenvalue decompositionのjulia実装がこちら．．

function reigen(A,r;p=10, twopass=false)
    n = size(A)[1]
    m = r + p
    Ω = randn(n, m)
    Y = A*Ω 
    if twopass
        Q = qr!(Y).Q * Matrix(I,n,m)
        T = Q' * A * Q
    else
         Q = qr(Y).Q * Matrix(I,n,m)
         T = (Q' * Y) * ( (Q' * Ω)\I)
    end
    λ, S = eigen(T)
    Λ = Diagonal(λ)
    U = Q*S'
    return U, λ
end

twopass=falseの方が速いが，より不正確です．twopass=trueの方がより正確です．pはover-sampling parameterで，この文献に従ってp=10としています．よい決め方知っている人がいたら教えてほしいです． ちょっと詰まったのがQR分解のところですね．qr(Y).Qはそのままでは正方行列なのでn×mの矩形行列として取り出したいときはqr(y).Q*Matrix(I,n,m)とします． 論文によって，QがQR分解のQだったり，特異値分解の主要な右特異ベクトルだったりしますが，これは，Johnson-Lindenstraussの補題で正当化できるって，強い人に教わりました．

なお，N Halko 論文によると，Y=A*Ωの積は，FFTを使うと，もっとはやく計算できるっぽいです．FFTとかあんまり分からないので，とりあえずこのままで進めます． 全体の計算量については，強い方に教えてもらいました．

精度について

さて，じゃぁ，計算時間の比較の前にeigenと同じ結果を得られるか実験してみます． 適当な30×30の対称行列Aの固有値をr=6個取り出してみます．

n = 30
r = 6
A = Symmetric(rand(n,n))

U, λ = reigen(A,r,twopass=false);
λ[end-r+1:end]
#  1.7535736114874219
#  2.2731278141844804
#  2.778424916550649
#  2.9919084354863856
#  8.475088453581433
# 14.989196447947041

U, λ = reigen(A,r,twopass=true);
λ[end-r+1:end]
#  1.2197855614132485
#  1.55690645146256
#  2.0387086285543066
#  2.111767689505717
#  2.3878222860931304
# 15.152842264762363

eigen(A).values[end-r+1:end]
# 1.7131238291530306
#  2.0291194695459063
#  2.2267398731934835
#  2.3956403884643858
#  2.6719356178806133
# 15.209793778533955

eigen(A)の出力結果と比べると，セルフ実装の結果が全然違うやん！！！ どこに誤りがあるのか分からなくて，StackOverflowで質問してしまったが，どこにも誤りはないです． というのも，上であげたどちらの論文でも説明されているように，REDは，固有値が減衰する行列に対して有用な手法のようです．知らなかった！！． ちなみに，固有値が減衰している場合もしていない場合も，RSVDによる行列の再構成に関しては，フロベニウスノルムの意味ではよく近似できるらしいが，特異値についてはあまりよく近似できていないっぽい．

If they decay fast, then the spectral norm error and the Frobenius norm error are similar,and the RSVD works well. If they decay slowly, then the RSVD performs fine when errors are measuredin the Frobenius norm, but not very well when the spectral norm is the one of interest. ([1607.01649] Randomized methods for matrix computations, P.G. Martinsson)

N halko論文の実験パートもよく読むと，固有値が減衰する手法についてRSVDを実験していますね．ということで，行列Aをいい感じに固有値が減衰する行列に書き直して再度実験すると，

n = 30
r = 6
A = Symmetric(rand(n,n))
B = zeros(n,n)
e, v = eigen(A)
for i in 1:n
    B += e[i] * v[:,i] * v[:,i]' * (i>n-r ? 1 : 1e-3)
end
A = Symmetric(B);

U, λ = reigen(A,r,twopass=false);
λ[end-r+1:end]
#  1.6972269946076446
#  2.0678582596060338
#  2.114020877070649
#  2.645103836439442
#  2.80298502359695
# 14.78544974888867

U, λ = reigen(A,r,twopass=true);
λ[end-r+1:end]
#  1.6883681640086117
#  2.0536447837869174
#  2.1103655428699253
#  2.6438950659892435
#  2.797125374675521
# 14.78174207245349

eigen(A).values[end-r+1:end]
#  1.6883689201446117
#  2.053646020688538
#  2.1103671910241286
#  2.643897161145393
#  2.7971278616415454
# 14.781742117602374

ちゃんとできてそう．

速度について

ということで，他の固有値分解の手法との計算速度を比較する実験をします．とりあえず，n=1500にして，固有値が100個目以降は小さくなるような行列を用意して， ターゲットランクrを徐々に大きくして計算時間の変化を観測しました．

using LinearAlgebra
using Arpack
using KrylovKit
using IterativeSolvers

n = 1500
A = Symmetric(rand(n,n))
B = zeros(n,n)
e, v = eigen(A)
for i in 1:n
    B += e[i] * v[:,i] * v[:,i]' * (i>n-100 ? 1 : 1e-3)
end
A = Symmetric(B);

p0 = 5
rs = [ p0*2^(n-1) for n = 1:7];

for r in rs
    result_lapack, runtime_lapack = @timed eigen(A, n-r+1:n);
    result_arpack, runtime_arpack = @timed eigs(A, nev=r, which=:LR);
    result_krylov, runtime_krylov = @timed eigsolve(A, r, :LM, krylovdim=r+10);
    result_lobpcg, runtime_lobpcg = @timed lobpcg(A, true, r);
    result_reigen1, runtime_reigen1 = @timed reigen(A, r, twopass=false);
    result_reigen2, runtime_reigen2 = @timed reigen(A, r, twopass=true);
    println(
       "r $r
        lapack  $runtime_lapack \t arpack $runtime_arpack \t krylov $runtime_krylov \t lobpcg $runtime_lobpcg
        reigen1 $runtime_reigen1 \t reigen2 $runtime_reigen2")
    println("******")
end

結果はこんな感じ．

r 5
        lapack  0.478838     arpack 0.0602754    krylov 0.0251191    lobpcg 0.1881867
        reigen1 0.0050092    reigen2 0.1158657
******
r 10
        lapack  0.5068946    arpack 0.0926673    krylov 0.0395786    lobpcg 0.3958909
        reigen1 0.0096395    reigen2 0.1367981
******
r 20
        lapack  0.5186003    arpack 0.0878039    krylov 0.0501138    lobpcg 0.4257776
        reigen1 0.0059776    reigen2 0.1848728
******
r 40
        lapack  0.5065934    arpack 0.1264882    krylov 0.0722774    lobpcg 0.8301309
        reigen1 0.0111938    reigen2 0.4549637
******
r 80
        lapack  0.794966     arpack 0.1150471    krylov 0.1482746    lobpcg 0.2263155
        reigen1 0.0181737    reigen2 0.7018366
******
r 160
        lapack  0.8155617    arpack 0.9227155    krylov 1.2557783    lobpcg 2.1098337
        reigen1 0.078549     reigen2 1.6126655
******
r 320
        lapack  0.8963038    arpack 3.7780156    krylov 3.7544088    lobpcg 6.2717872
        reigen1 1.1731766    reigen2 3.3344803

r<<nだとREDの圧勝ですね．めでたしめでたし． さて，私の研究では，分解したい行列Aの固有値は別に減衰していなさそうなので，REDは使えなそう．代わりにNystrom手法による固有値分解の手法とかでも実装しようかなと思います．

Julia言語 実対称行列の固有値を大きい順にr個抽出する．

Julia標準搭載のeigenを使うと，N×N行列Aの固有値分解がO(N^3)になっちゃうと思ってたからArpack使ってたんだけど，Arpackはそもそも疎行列向きっぽい． どうやら入力行列Aが実対称行列ならば

val, vec = eigen(A, 1:r)

で小さい順に固有値をr個とりだせるみたい．大きい順に取り出したい場合は

val, vec = eigen(A, N-r+1:N)

でいけるっぽい．ついでにArpackでは，

val, vec=eigs(A, nev=r, which=:LR);

で大きい順にrの固有値を計算する．

ということで，N=3500にして，rをうごかしてみて，計算時間がどんな感じになるかちょっと試してみた．

using LinearAlgebra
using Arpack

N = 3500
A = Symmetric(rand(N,N));

@time _,_ = eigen(A, N-5+1:N);
@time _,_ = eigen(A, N-10+1:N);
@time _,_ = eigen(A, N-20+1:N);
@time _,_ = eigen(A, N-40+1:N);
@time _,_ = eigen(A, N-80+1:N);
@time _,_ = eigen(A, N-160+1:N);
@time _,_ = eigen(A, N-320+1:N);
@time _,_ = eigen(A, N-640+1:N);
@time _,_ = eigen(A, N-1280+1:N);
@time _,_ = eigen(A, N-2560+1:N);

すると，こんな感じ．やっぱり，rを大きくすると徐々に計算量が増えていく．でも理論通りなら，計算量はO(N² r)なので，倍々で計算時間が増えていくと思うんだけどなぁ，おかしいな．

  2.236443 seconds (21 allocations: 138.497 MiB, 1.99% gc time)
  2.287247 seconds (21 allocations: 138.612 MiB)
  2.423010 seconds (21 allocations: 138.841 MiB)
  2.569268 seconds (21 allocations: 139.299 MiB)
  2.423246 seconds (21 allocations: 140.215 MiB)
  2.690961 seconds (21 allocations: 142.046 MiB, 0.28% gc time)
  3.180988 seconds (21 allocations: 145.710 MiB)
  4.213903 seconds (21 allocations: 153.036 MiB, 0.15% gc time)
  9.066842 seconds (21 allocations: 167.690 MiB)
 26.820623 seconds (20 allocations: 196.996 MiB, 0.03% gc time)

ちなみに

@time _,_ = eigen(A);

だと， 5.097789 seconds (19 allocations: 207.070 MiB, 0.12% gc time) でした．つまり，rがある程度大きい場合は，全部求めちゃった方が速かったりするってことですね．

次に，Arpackの結果はこんな感じ．

@time _,_ = eigs(A,nev=5, which=:LR);
@time _,_ = eigs(A,nev=10, which=:LR);
@time _,_ = eigs(A,nev=20, which=:LR);
@time _,_ = eigs(A,nev=40, which=:LR);
@time _,_ = eigs(A,nev=80, which=:LR);
@time _,_ = eigs(A,nev=160, which=:LR);
@time _,_ = eigs(A,nev=320, which=:LR);
@time _,_ = eigs(A,nev=640, which=:LR);
@time _,_ = eigs(A,nev=1280, which=:LR);
@time _,_ = eigs(A,nev=2560, which=:LR);

  0.548882 seconds (2.33 k allocations: 775.391 KiB)
  1.044984 seconds (4.04 k allocations: 983.375 KiB)
  0.806430 seconds (3.35 k allocations: 1.632 MiB)
  1.210613 seconds (4.00 k allocations: 3.070 MiB)
  1.765235 seconds (5.56 k allocations: 6.031 MiB)
  3.085386 seconds (8.10 k allocations: 12.222 MiB)
  8.542838 seconds (12.02 k allocations: 25.732 MiB)
 24.201966 seconds (16.58 k allocations: 57.314 MiB)
 101.972415 seconds (23.18 k allocations: 139.133 MiB)
  97.998941 seconds (18.07 k allocations: 196.968 MiB, 0.03% gc time)

いまAは密行列ですが，rが小さいときは，疎行列用のArpackの方が速いという結論ぽいですね．ちなみにeigen(A)とすると，デフォルトで6個だけ固有値を抽出する仕様っぽいです．

ちなみに，KrylovKit.jlというライブラリもあって，こちらは

using KrylovKit
@time _,_,_ = eigsolve(A, 5, :LR);
@time _,_,_ = eigsolve(A, 10, :LR);
@time _,_,_ = eigsolve(A, 20, :LR);

 0.322175 seconds (2.65 k allocations: 15.157 MiB)
 0.371782 seconds (3.62 k allocations: 20.360 MiB)
 0.899626 seconds (12.60 k allocations: 71.841 MiB)

という風にeigsolve(A, r, :LR)でr個の固有値を取り出します．rを大きくすると，krylovdimが小さいとエラーがでるのですが，krylovdimをどれくらい大きくするのが最適なのかよくわからなかったです．

Julia言語 データをTrainとValidにsplitする．

Juliaで，n×m行列Xに，n次元の特徴量がm個格納されていて，これを，trainとvalidに分けたい． 要は，列をランダムにk個選んだn×k行列とn×(m-k)行列をつくりたい．

手っ取り早い方法がわかんないんだけど，いまのところ，暫定のマイベストプラクティスはこんな感じ．

using InvertedIndices
using StatsBase

function data_split(X, ratio)
    n, m = size(X)
    k = Int(floor(ratio*m))
    idxs = sample(1:m, k, replace=false, ordered=true)
    X_train = X[:, idxs]
    X_valid = X[:, Not(idxs)]
    return X_train, X_valid
end

教師データYも一緒にsplitしたいなら，こんな感じ．

using InvertedIndices
using StatsBase

function data_split(X, Y, ratio)
    n, m = size(X)
    k = Int(floor(ratio*m))
    idxs = sample(1:m, k, replace=false, ordered=true)
    X_train = X[:, idxs]
    Y_train = Y[:, idxs]
    X_valid = X[:, Not(idxs)]
    Y_valid = Y[:, Not(idxs)]
    return X_train, Y_train, X_valid, Y_valid
end

もっと良い方法があったら知りたい．

Julia言語でカーネル回帰

カーネル回帰のjulia実装が研究で必要になったから１０分で実装した．

とりあえず，前に書いたカーネルを呼ぶ式．

function rbf_inners(X,Y,c)
    rbf_inner(x,y,c) = exp( -c*norm(x-y)^2 )
    d, n = size(X)
    d, m = size(Y)
    K = zeros(m,n)
    for j = 1:n
        for i = 1:m
           K[i,j] = rbf_inner(Y[:,i],X[:,j],c)
        end
    end
    return K
end

function poly_inners(X,Y,c,deg)
    poly_inner(x,y,c,deg) = (c+x'*y)^deg
    d, n = size(X)
    d, m = size(Y)
    K = zeros(m,n)
    for j = 1:n
        for i = 1:m
           K[i,j] = poly_inner(Y[:,i],X[:,j],c,deg)
        end
    end
    return K
end

function lapras_inners(X,Y,c)
    lapras_inner(x,y,c) = exp(-c*sum( abs.(x .- y)))
    d, n = size(X)
    d, m = size(Y)
    K = zeros(m,n)
    for j = 1:n
        for i = 1:m
           K[i,j] = lapras_inner(Y[:,i],X[:,j],c)
        end
    end
    return K
end

function Kernel(X, Y; c=1.0, deg=3.0, kernel="rbf")
    if kernel == "rbf"
        K = rbf_inners(X,Y,c)
    elseif kernel == "poly"
        K = poly_inners(X,Y,c,deg)
    elseif kernel == "linear"
        K = poly_inners(X,Y,0,1)
    elseif kernel == "lapras"
        K = lapras_inners(X,Y,c)
    else
        error("Kernel $kernel is not defined")
    end
    return K
end

カーネル回帰の本体はこんな感じ．正則化のパラメータがlambda．これが大きいと正則化が強くなる． cとdegがカーネルのパラメータ．cはRBFカーネルの定数で，dは多項式カーネルの次数． Lは訓練データの数．逆行列を求めるなら，\の方がいいってどこかで誰かが言っていた気がするので，invは使ってない． 機械学習民なら，カーネル回帰の式は特に何の説明もいらないと思うがわかりやすい解説はQiitaにあると思う．（これとか）

function Kreg(X_train, Y_train, X_test; lambda=1, c=1, deg=3.0, kernel="rbr")
    K = Kernel(X_train, X_train)
    L = size(K)[1]
    alpha = (K + lambda * Matrix(1.0I, L, L)) \ Y_train'
    k = Kernel(X_train, X_test)
    Y_pred = k * alpha
    return Y_pred
end

適当にデータを作ってプロットしてみる．

using LinearAlgebra
using Distributions
function make_1d_data(L)
    f(x) = x/2 + sin(x)*cos(x)+1
    dist = Uniform(-6, +6)
    x = rand(dist,(1, L))

    eps = 0.1
    noise_dist = Uniform(-eps, +eps)

    y = zeros(1, L)
    for i = 1:L
        noise = rand(noise_dist)
        y[1,i] = f(x[i]) + noise
    end
    return x, y
end

function main()
    L = 150
    N = 50
    X_train, Y_train = make_1d_data(L)
    X_test, Y_test = make_1d_data(N)
    Y_pred = Kreg(X_train, Y_train, X_test)
    L2error = norm(Y_test[:] .- Y_pred[:])^2 / size(Y_test)[2]
    println("Ferror:$L2error")
    plots(X_train, Y_train, X_test, Y_pred)
end

あっているように思います．ハイパラcとlambdaを調節したらもっと汎化すると思う． ちなみにプロットに使ったのはこちらのコード

using Plots
using Plots.PlotMeasures
function plots(x, y, x_test, y_pred)
    fnt = font(20, "Helvetica")
    plt = plot(legend=:topleft, legendfont=fnt,
                bottom_margin = 10mm,
                left_margin = 10mm,
                size=(900,500),
                xlabel ="x", ylabel="y",
                xguidefont = fnt,
                yguidefont = fnt,
                xtickfont = fnt,
                ytickfont =fnt
              )
    plot!(plt, x[:], y[:], label="train", color="blue", st=:scatter)
    plot!(plt, x_test[:], y_pred[:], label="test", color="red", st=:scatter)
    png("tmp.png")
    println("saved")
end

追記

このカーネル回帰の実装は，重回帰にしてもちゃんと動く． ためしに，下の人工データを使ってみるとちゃんとコスト関数が小さくなってる．

function make_2d_data(L)
    #f(x1,x2) = sin(x1+x2^2)
    f(x1,x2) = (x1^2+x2^2)*exp(1-x1^2-x2^2)
    dist = Uniform(-2, +2)
    x = rand(dist,(2, L))

    y = zeros(1, L)
    for i = 1:L
        y[1,i] = f(x[1,i],x[2,i])
    end
    return x, y
end

Ferror:0.008802513523177683

Julia言語 LightGBMで多クラス分類

（解決済み）

非深層学習系の提案手法を実装するならJulia一択だと思うのですが，既存の比較手法をすぐに試せないとpython民に戻りたくなりますね．．LightGBMで多クラス分類をしようとして，ちょっと詰まったので，二度と再び本件で躓かないためのマイメモです．

LightGBMはここにあるgithubのものを落としてきました．

github.com

まぁ普通に]add LightGBMでインストールできるのですが，提供されているデモは二値分類用．多クラス分類のためにはいろいろ変えないといけない．

using LightGBM
using DelimitedFiles

function LGBM(x_train, y_train, x_test, y_test, num_class)
    y_test = convert(Vector, Int.(y_test))
    y_train = convert(Vector, Int.(y_train))

    y_test .-= 1 # labelは0からスタートっぽいので，
    y_train .-= 1

    # Load LightGBM's binary classification example.
    # Create an estimator with the desired parameters—leave other parameters at the default values.
    estimator = LGBMClassification(
        objective = "multiclass",
        num_iterations = 100,
        learning_rate = .1,
        early_stopping_round = 5,
        feature_fraction = .8,
        bagging_fraction = .9,
        bagging_freq = 1,
        num_leaves = 1000,
        num_class = num_class,
        metric = ["auc_mu", "multi_logloss"]
    )

    LightGBM.fit!(estimator, x_train, y_train)
    predicted = LightGBM.predict(estimator, x_test)
end

はまりポイント① なぜかfit!やpredictが定義されてない，って怒られるので，LightGBM.fit!とかLightGBM.predictとする．

はまりポイント② metricのaucはニクラス分類用．多クラス分類をするときは，auc_muとしよう．同様に，コスト関数にはmulti_loglossを使いましょう．

はまりポイント③ ラべルyの各要素はIntにしておこう．ラベルの始点は0なので，y_test, y_trainが1からはじまるようにロードしているなら，.-1しましょう．

とりあえずこれで動いた！！ predictedは[テストサンプル数，クラス数]の二次元配列になっていて，各サンプルがどのクラスにいるのかの確率をあらわしてるっぽい．

Julia言語 Kernel PCA実装

カーネルPCAを実装した．参考にしたのは 赤穂先生のカーネル本とこの修士論文．

「RKHS内でのΦ(x)の平均を簡単のためにゼロにする」という巷の説明をあちこちで見かけるが，無限次元空間でΦ(x)の平均をゼロにする方法は自明じゃない．この修論のappendixがとても参考になった．

https://amzn.to/3zFF2wV

まず，カーネル行列を計算するための関数を書く． （この辺，後々，testデータとtrainデータ間のカーネルも必要になったりするので，地味に書き方が難しい．カーネル法のライブラリってこの辺，どうしてるんだろう）

using LinearAlgebra

function rbf_inners(X,Y,c)
    rbf_inner(x,y,c) = exp( -c*norm(x-y)^2 )
    d, n = size(X)
    d, m = size(Y)
    K = zeros(m,n)
    for j = 1:n
        for i = 1:m
           K[i,j] = rbf_inner(Y[:,i],X[:,j],c) 
        end
    end
    return K
end

function poly_inners(X,Y,c,deg)
    poly_inner(x,y,c,deg) = (c+x'*y)^deg
    d, n = size(X)
    d, m = size(Y)
    K = zeros(m,n)
    for j = 1:n
        for i = 1:m
           K[i,j] = poly_inner(Y[:,i],X[:,j],c,deg) 
        end
    end
    return K
end

function lapras_inners(X,Y,c)
    lapras_inner(x,y,c) = exp(-c*sum( abs.(x .- y)))
    d, n = size(X)
    d, m = size(Y)
    K = zeros(m,n)
    for j = 1:n
        for i = 1:m
           K[i,j] = lapras_inner(Y[:,i],X[:,j],c) 
        end
    end
    return K 
end

function Kernel(X, Y; c=2.0, deg=3.0, kernel="rbf")
    if kernel == "rbf"
        K = rbf_inners(X,Y,c) 
    elseif kernel == "poly" 
        K = poly_inners(X,Y,c,deg)
    elseif kernel == "linear"
        K = poly_inners(X,Y,0,1)
    elseif kernel == "lapras"
        K = lapras_inners(X,Y,c)
    else
        error("Kernel $kernel is not defined")
    end
    return K
end

次に，カーネルPCAの本体を書いていく．KからKbarを得る更新式は，‪Bernhard Schölkopfさんの修論のappendixに書いてある．良い性質のカーネルを選ばないと，Φ空間での特徴量ベクトルの平均が0になっていることはない(と思う)ので，centered_processedは基本的にはfalseにするのが良いと思う．

rankがターゲットランクで，無限次元の特徴量空間につくる有限次元の線形部分空間の次元に対応する．訓練に使ったデータを，Φ空間のrank次元部分空間に射影した先のrank次元の座標がprojected_lowdim_pointsになる．射影した点に興味はなくて，基底だけほしい場合は，projection=falseにする．

function Kernel_PCA(X ; rank=2, kernel="rbf", c=2.0, deg=3.0, centered_processed=false, projection=true)
    K = Kernel(X, X, c=c, deg=deg, kernel=kernel)
    n, n = size(K)
    if centered_processed == true 
        Kbar = K
    else
        one = ones(n,n)
        Kbar = K - 1/n * one * K - 1/n * K * one + 1/(n^2) * one * K * one
    end
    Kbar = Symmetric(Kbar)
    
    f = eigen(Kbar, sortby=-)
    lambda = f.values ./ n
    alpha = f.vectors
    
    projected_lowdim_points = zeros(rank,n)
    if projection
        projected_lowdim_points = n * Diagonal(lambda[1:rank]) * alpha[:,1:rank]' 
    end
    return K, lambda, alpha, projected_lowdim_points
end

途中のSymmetricはなくてもよいはずだが，いれないと，微妙に固有値が複素数になることがあるので，念のため．基本的にはこれで実装完了で，ちょっと実験してみる．とりあえず適当に線形分離できない2次元データを作ってプロットしてみる．

using Plots
X = zeros(2,200)
for j = 1:100
    X[:,j] = [cos(2*pi*j/100.0)+0.1*rand() sin(2*pi*j/100.0)+0.1*rand()]
end
for j = 101:200
    X[:,j] = [0.25*cos(2*pi*j/100.0)+0.1*rand() 0.25*sin(2*pi*j/100.0)+0.1*rand()]
end

plt=scatter(X[1,1:100],X[2,1:100], c="blue", legend=false)
scatter!(plt,X[1,101:200], X[2,101:200], c="red")

もちろん，こんなのは，古典的なPCAじゃ手に負えない．これをKernel-PCAして，Φ空間での2次元空間にプロットしてみる．

K, lambda, alpha, lowdims=Kernel_PCA(X, rank=2, kernel="rbf", c=2.0);
plt=scatter(lowdims[1,1:100],lowdims[2,1:100], c="blue", legend=false)
scatter!(plt,lowdims[1,101:200], lowdims[2,101:200], c="red")

素晴らしい！！！ きちんと線形分離できとる！！

データXに加えて，データXtestが観測された場合を考える． 素朴に考えれば，XとXtestをconcateして，同じことを繰り返せばよいが，それは手間である．固有値を求めるのは計算コストが高いので同じことを繰り返したくない．．．そこで，やはり彼の修論のappendixを参照にして，fitさせる関数をつくる．

function fit_Kernel_PCA(K, Xtrain, Xtest, alpha; rank, kernel="rbf", c=2.0, deg=3.0, centered_processed=false)
    Ktest = Kernel(Xtrain, Xtest, c=c, deg=deg, kernel=kernel)
    m, l = size(Ktest)
    if centered_processed 
        Ktest_bar = Ktest
    else
        one_ml = ones(m,l)
        one_ll = ones(l,l)
        Ktest_bar = Ktest - 1/l * Ktest*one_ll - 1/l * one_ml * K + 1/(l^2) * one_ml * K * one_ll
    end
    projected_lowdim_points = (Ktest_bar * alpha[:,1:rank])' 
    return projected_lowdim_points
end

これを使うと，固有値計算を繰り返す必要はない（テストサンプルとトレーニングサンプル間のカーネルを計算する必要はあるが，これは，固有値計算に比べればずっと軽い）．実験してみよう．×印が，新たに観測されたデータで，〇印が，Φでrank次元線形部分空間をつくるのに使ったデータ．

Xtest = zeros(2,100)
for j = 1:50
    Xtest[:,j] = [cos(2*pi*j/50.0)+0.1*rand() sin(2*pi*j/50.0)+0.1*rand()]
end
for j = 51:100
    Xtest[:,j] = [0.25*cos(2*pi*j/50.0)+0.1*rand() 0.25*sin(2*pi*j/50.0)+0.1*rand()]
end

plt=scatter(X[1,1:100],X[2,1:100], c="blue", legend=false)
scatter!(plt,X[1,101:200], X[2,101:200], c="red")
scatter!(plt,Xtest[1,1:50], Xtest[2,1:50], c="blue", markershape=:cross)
scatter!(plt,Xtest[1,51:100], Xtest[2,51:100], c="red", markershape=:cross)

lowdims_test =fit_Kernel_PCA(K, X, Xtest, alpha; rank=2, kernel="rbf", c=2.0, deg=3.0, centered_processed=false)
plt=scatter(lowdims[1,1:100],lowdims[2,1:100], c="blue", legend=false)
scatter!(plt,lowdims[1,101:200], lowdims[2,101:200], c="red")
scatter!(plt,lowdims_test[1,1:50],lowdims_test[2,1:50], c="blue", markershape=:cross)
scatter!(plt,lowdims_test[1,51:100],lowdims_test[2,51:100], c="red", markershape=:cross)

正しく実装できてる．

Python αダイバージェンスを小さくするNMF

KLダイバージェンスを最小化するNMFはよく知られていますし，sklearnで'beta_loss=1'にすればすぐに利用できます． 今回はKLダイバージェンスの拡張であるαダイバージェンスを小さくするNMFを実装しました．

この論文の更新式をさくっと実装しました． www.sciencedirect.com

収束条件とかはsklearnの普通のNMFに揃えようとしたけど，ちょっと違うかも． （もしかしたら，sklearnの方はコスト関数のルートをとっているかもしれない）

def D_KL(P,Q):
    epsilon = 1.0e-10
    P = P+epsilon
    Q = Q+epsilon
    return np.sum(P*np.log(P/Q)) - np.sum(P) + np.sum(Q)

def alpha_KL(X,AS,alpha):
    if alpha == 1.0:
        return D_KL(X, AS)
    if alpha == 0.0:
        return D_KL(AS, X)
    
    alpha_kl = 0.0
    n, m = np.shape(X)
    term = 0.0
    for i in range(n):
        for j in range(m):
            term += alpha * X[i,j] + (1-alpha) * AS[i,j] - math.pow(X[i,j],alpha) * math.pow((AS[i,j]),1-alpha)
    
    alpha_kl = term / (alpha*(1-alpha))
    
    return alpha_kl
    
def nmf_alpha(X, r, alpha, cnt_max=100, tol=1.0e-3, verbose=False):
    n, m = np.shape(X)
    S = np.random.rand(r, m) 
    A = np.random.rand(n, r)
    
    if alpha == 0.0:
        alpha += 1.0e-5
    
    error_at_init = alpha_KL(X, np.dot(A,S),alpha)
    previous_error = error_at_init
    for cnt in range(cnt_max):
        # update S
        for i in range(r):
            for j in range(m):
                sumA = np.sum(A, axis=0)
                term = 0.0
                AS = np.dot(A,S)
                for k in range(n):
                    term += A[k,i] * math.pow(( X[k,j] / AS[k,j] ),alpha)
                S[i,j] *= pow( term / sumA[i], 1.0/alpha )

        # update A
        for i in range(n):
            for j in range(r):
                sumS = np.sum(S, axis=1)
                term = 0.0
                AS = np.dot(A,S)
                for k in range(m):
                    term += S[j,k] * math.pow(( X[i,k] / AS[i,k] ),alpha)
                A[i,j] *= pow( term / sumS[j], 1.0/alpha )
                
        if tol>0 and cnt % 10 == 0:
            error = alpha_KL(X, np.dot(A,S), alpha)
            if verbose:
                print(f"iter:{cnt} cost:{error}")
            if (previous_error - error) / error_at_init < tol:
                break
            previous_error = error
    return A, S

とりあえず，こんな感じで実行． 分解後のターゲットランクは2にして，αは1.0にしてみる．

X = np.random.rand(30,30)
A,S = nmf_alpha(X, 2, 1.0, verbose=True, cnt_max=200)

実行結果はこんな感じ．

iter:0 cost:88.0270095179497
iter:10 cost:77.71017624081105
iter:20 cost:75.81128050686453
iter:30 cost:73.81279529499801
iter:40 cost:72.67141167025983
iter:50 cost:72.32208831953886
iter:60 cost:72.20333908949885

ちゃんとロスがさがっているので正しいっぽい． α=1の極限では，普通のKL-NMFと一致するので，sklearnの結果と一致するか検証． beta_loss=1はコスト関数をKLにするということ．

from sklearn.decomposition import NMF
sklearn_NMF = NMF(n_components=2, init='random', random_state=0, max_iter=200, beta_loss=1, solver='mu')
W = sklearn_NMF.fit_transform(X)
H = sklearn_NMF.components_
R = np.dot(W,H)
cost = alpha_KL(X, R,1.0)
print(cost)

すると

72.21912349939544

まぁ，ロスが同じくらいなのであってそう．

Python 確率単体上で凸関数を最小化(SciPy)

多変数関数f(x)を確率単体上のxで最適化したい．（最近，情報幾何をやっているので，こういうタスクばかりやってる気がする） つまり，xは 0 < x[i] < 1 かつ sum(x) = 1 の範囲でf(x)の最小化をしたい． ちまちま自分でソルバー書いてたけど，scipy使う方がよさそう． ドキュメント見た感じ，これでいける．

from scipy import optimize

def f(x):
    return (x[0]-1)**2 + (x[1]-1)**2 + (x[2]-1)**2

cons = [ 
        {'type':'eq', 'fun': lambda x: 1-sum(x)},
        {'type':'ineq', 'fun': lambda x: -x+1},
        {'type':'ineq', 'fun': lambda x: x}
       ]

optimize.minimize(f, np.array([-3,-3,-3]), constraints=cons)

Julia言語 benchmark で計測した時間を変数に代入したい．

BenchmarkToolsって便利ですよね．例のように(10×10)のランダム行列の固有値を計算するのに要した時間を計測したとします．

ここから，最小時間37.748とか，中央値47.454とかを取り出すにはどうしたらよいでしょうか，というのがいつも分からなくなるのでメモ． おそらく，最適解は以下のように，time(median(t))とするのが良いと思う．これは強い人に教わった．time() で取得できるのがナノ秒単位なので，time(median(t))/1000で，μ秒になります．

ちなみに，time(std(t))のためには，using Statistics をしておく必要があります．（はまった．．．）

また，こういうときの対処法として，このツイートがとても参考になります．

Julia言語 CartesianIndexからvecに変換する方法

しょっちゅう分からなくなってはTwitterで強い人に教えてもらうので，ここにメモ．

CartesianIndexは.Iで添え字を取り出せる．これを'collect'すると所望のvecになる.

ただ，これだとちょっと気持ち悪いと思う人もいるようなので，

Tuple(CartesianIndex(2,3,4))

でタプルにして，これをcollectする

collect(Tuple(CartesianIndex(2,3,4)))

のが良いと強い人に教わりました．

タプルなど多くのものを配列にするにはcollectが便利のようだ！

Julia言語 補助行列が2つある場合のNon-negative Multiple Matrix Factorization (NM2F)の実装

前回，補助行列が一つだけのNMMFの実装を公開した． 今回は同様に，これらの論文にのっとり，補助行列が２つある場合のNMMFを実装した． 行列X,Y,Zを，WH, AH,WBに同時に分解するというタスクである．コスト関数はKL情報量．（LSの場合の更新式ってどこかにない？？）

www.ijcai.org

ci.nii.ac.jp

原著論文に忠実に，一様離散分布で初期化する実装は以下．

function cost(X,Y,Z,W,H,A,B,alpha,beta)
    d_kl(P,Q) = sum( P .* log.( P ./ Q ) ) - sum(P) + sum(Q)
    total_cost = d_kl(X, W*H) + alpha*d_kl(Y, A*H) + beta*d_kl(Z, W*B)
    return total_cost
end

"""
Non-negative Multiple Matrix Factorization (NM2F)
See [original paper 1](https://www.ijcai.org/proceedings/2017/407),
[original paper 2](https://ci.nii.ac.jp/naid/110009691643)
# Aruguments
- 'X' target matrix. Factorized to the bases W and the coefficients H.
- 'Y' auxiliary feature. Factorized to A*H
- 'Z' auxiliary data. Factorized to W*B
- 'K' target rank, int.
- 'alpha' scale of auxiliary feature
- 'beta'  scale of auxiliary data
"""
function NM2F(X, Y, Z, K; alpha=1.0, beta=1.0, max_iter = 200, tol = 1.0E-4, verbose = true)
    W = rand(I, K)
    H = rand(K, J)
    A = rand(N, K)
    B = rand(K, M)
    Xhat = W*H
    Yhat = A*H
    Zhat = W*B
    
    error_at_init = cost(X,Y,Z,W,H,A,B,alpha,beta)
    previous_error = error_at_init
    for iter = 1 : max_iter
        W .*= ( (X ./ Xhat) * H' .+ beta * (Z ./ Zhat) * B') ./ ( ones(I,J) * H' .+ beta * ones(I, M) * B'  )
        Xhat = W*H
        Zhat = W*B
        
        H .*= ( W' * (X ./ Xhat) .+ alpha * A'*(Y ./ Yhat) ) ./ ( W' * ones(I,J) .+ alpha * A' * ones(N, J) )
        Xhat = W*H
        Yhat = A*H
        
        A .*= ( (Y ./ Yhat) * H'  ) ./ ( ones(N,J) * H'  )
        Yhat = A*H
        
        B .*= (  W' * (Z ./ Zhat) ) ./ ( W' * ones(I, M) )
        Zhat = W*B
        
        error = cost(X,Y,Z,W,H,A,B,alpha,beta)
        if verbose && iter % 20 == 1
            println("iter:$iter cost:$error")
        end
        if (previous_error - error) / error_at_init < tol
            break
        end
        previous_error = error
    end
    return W, H, A, B
end

I = 6
J = 5
M = 4
N = 3

X = rand(I, J)
Y = rand(N, J)
Z = rand(I, M)
W,H,A,B = NM2F(X,Y,Z,2);

Kはターゲットランク． このプログラムを実行すると，コスト関数がどんどん下がっていく様子が分かる．

iter:1 cost:6.1250592119187575
iter:21 cost:4.615593687122641
iter:41 cost:4.465480179778419
iter:61 cost:4.150666236462589
iter:81 cost:3.854919025274069

Julia言語 Non-negative Multiple Matrix Factorization (NMMF)の実装

複合データ分析技術についてこのpdfで勉強してたら，NMMFというものを知った．

ちゃんとした議論は，IJCAIでしている． www.ijcai.org

人工知能学会の予稿や論文，絵本のパターン抽出や，和文論文が参考になる．

私も研究で使うことになったのでささっと実装した．

とりあえず，いつも使っているコスト関数

using LinearAlgebra

function get_dkl(P, Q ; eps_thred = 1.0e-10)
    tensor_shape = size(P)
    dkl = 0.0
    for idx in CartesianIndices(tensor_shape)
        if P[idx] < eps_thred
            continue
        end
        if Q[idx] < eps_thred
            continue
        end
        dkl += P[idx] * log( P[idx] / Q[idx] )
    end
    return dkl - sum(P) + sum(Q)
end

function cost_function(cost, X, Xr)
    if cost == "LS"
        return norm(X-Xr)

    elseif cost == "KL"
        dkl = get_dkl(X, Xr)
        return dkl
    end
end

としておく． 本体は以下．これで実行すれば，ロスが下がっていく様子がわかる． 初期化は超適当に連続一様分布から行っている．

"""
Non-Negative Multiple Matrix Factorization
proposed by Koh Takeuchi, Katsuhiko Ishiguro, Akisato Kimura and Hiroshi Sawada
See [original paper](https://www.ijcai.org/Proceedings/13/Papers/254.pdf)
You can see the implementation in [this Japanese documents]
(http://www.kecl.ntt.co.jp/icl/ls/members/tatsushi/PDF/IEICE_vol99_no6_543-550.pdf)

# Aruguments
- `X` input I \times J matrix
- `Y` input I \times K matrix
- `R` target rank
"""
function NMMF_KL(X, Y, R; max_iter=200, tol = 1.0E-3, verbose = false)
    (I, J) = size(X)
    (I, K) = size(Y)
    A = rand(I, R)
    B = rand(J, R)
    C = rand(K, R)

    Xhat = A * B'
    Yhat = A * C'

    cost_at_init  = cost_function("KL", X, Xhat) + cost_function("KL", Y, Yhat)
    previous_cost = cost_at_init
    for iter = 1 : max_iter
        A .= A .* ( ( (X ./ Xhat ) * B + (Y ./ Yhat) * C) ./ ( ones(I,J) * B + ones(I,K) * C ) )

        Xhat = A * B'
        Yhat = A * C'

        B .= B .* (( (X ./ Xhat)' * A ) ./ (ones(J,I)*A))
        Xhat = A * B'

        C .= C .* (( (Y ./ Yhat)' * A ) ./ (ones(K,I)*A))
        Yhat = A * C'

        if tol > 0 && (iter % 10 == 0)
            cost = cost_function("LS", X, Xhat) + cost_function("LS", Y, Yhat)
            if verbose
                println("step:$iter loss:$cost")
            end
            if (previous_cost - cost) / cost_at_init < tol
                break
            end
            previous_cost = cost
        end
    end
    return A, B, C
end

function NMMF_LS(X, Y, R; max_iter=200, tol = 1.0E-3, verbose = false)
    (I, J) = size(X)
    (I, K) = size(Y)
    A = rand(I, R)
    B = rand(J, R)
    C = rand(K, R)

    Xhat = A * B'
    Yhat = A * C'

    cost_at_init  = cost_function("LS", X, Xhat) + cost_function("LS", Y, Yhat)
    previous_cost = cost_at_init
    for iter = 1 : max_iter
        A .= A .* ( (X * B + Y * C) ./ ( Xhat * B + Yhat * C ))
        Xhat = A * B'
        Yhat = A * C'

        B' .= B' .* ( (A' * X) ./ ( A' * Xhat ) )
        Xhat = A * B'

        C' .= C' .* ( (A' * Y) ./ ( A' * Yhat ) )
        Yhat = A * C'

        if tol > 0 && (iter % 10 == 0)
            cost = cost_function("LS", X, Xhat) + cost_function("LS", Y, Yhat)
            if verbose
                println("step:$iter loss:$cost")
            end
            if (previous_cost - cost) / cost_at_init < tol
                break
            end
            previous_cost = cost
        end
    end
    return A, B, C
end

X = rand(20,10)
Y = rand(20,8)
A,B,C = NMMF_LS(X, Y, 2, verbose=true)
A,B,C = NMMF_KL(X, Y, 4, verbose=true)

Julia言語 Non-Negative Tucker Decomposition の実装 [完全版] LSエラー, KLエラー

昨日の記事の完全版．Nonnegative Tucker Decomposition の Julia 実装． 重要なアルゴリズムと思うけど，Julia版は探しても見つからなかったので，自分で書いた．

ieeexplore.ieee.org

Python版はすでに実装があるみたい．これはSVDで初期化しているっぽいね．

github.com

論文中にいくつかトラップがある．（追記）

・式(17)は間違い． A^n-1 ... A¹ \Otimes A^N ... Aⁿ⁺¹が正しい．（まぁ，これは実装には関係ないけど．）

・Table4の S←Rand(R,m)はS←Rand(R,l)が正しいと思われる．

・Table4, Table5に無定義で出てくる行列の肩の＋は転置Tの誤植と思われる．（じゃないと行列の積が演算できない）

・更新は，A¹→A²→…→A^N→S→A¹...の順でやる．A¹→S→A²→S→A³→....なのか迷ったけど，論文にはIn updating A, the rest of parameters denoted by Sn are fixedとあるので，前者で正しいと思う．

・式(27)のεは入力Xと同じサイズのテンソル．

あと，式(22)を使うと計算効率がよくなるらしいが，これ，両辺，行列のサイズ違くないか？ 何度か計算したがつじつまが合いそうになかったで，とりあえず式(22)を使わない実装をした．（僕が勘違いしているようだったらコメントで教えてほしい）

初期化アルゴリズムは，論文に倣った”NMF"と"Random"の二つを用意した． 適当にrand(48,48,24)でテンソルつくって，(24,24,30)にランクを落としてみたら，NMFで初期化するとコスト関数は初めからかなり下がりきっているようだった．Figure1.を見ると，Iterations in initialization step are also added in these error traceとある．Figure1.のstep 1-50 くらいでかなりコスト関数が落ちているように思うので，初期化でかなり収束に近づいているっぽい？

using LinearAlgebra
using TensorToolbox

function cost_function(cost, X, Xr)
    if cost == "LS"
        return norm(X-Xr)
    elseif cost == "KL"
        return sum( X .* log.( X ./ Xr ) ) - sum(X) + sum(Xr)
    end
end

function reconst(S, A)
    N = ndims(S)
    Xr = S
    for n=1:N
        Xr = ttm(Xr, A[n], n)
    end
    return Xr
end

function init_random(tensor_size, reqrank)
    A = []
    N = length(reqrank)
    for n = 1:N
        An = rand(tensor_size[n], reqrank[n])
        push!(A, An)
    end
    S = rand( reqrank... )
    return S, A
end

"""
NMF Initializer. See Table 4 in the paper.
"""
function init_for_NMF(X, R; eps=1.0e-10, pre_iter_max_inloop=40)
    (m, l) = size(X)
    A = rand(m, R)
    S = rand(R, l)

    for iter = 1:pre_iter_max_inloop
        A .= max.(X*S', eps)
        A .= A .* ((X*S') ./ (A*S*S'))
        S .= max.(A'*X, eps)
        S .= S .* ( (A'*X) ./ (A'*A*S) )
    end
    return S, A
end

"""
Initializer based on NMF. See Table 5 in the paper.
"""
function init_based_NMF(X, reqrank; eps=1.0e-10, pre_iter_max=40, pre_iter_max_inloop=40)
    N = ndims(X)
    A = []
    for n=1:N
        _, An = init_for_NMF(tenmat(X, n), reqrank[n], pre_iter_max_inloop=pre_iter_max_inloop)
        push!(A, An)
    end

    S = rand( reqrank... )
    for iter = 1:pre_iter_max
        for n=1:N
            SAn = S
            for m in [1:n-1; n+1:N]
                SAn = ttm(SAn, A[m], m)
            end
            SAn = tenmat(SAn, n)
            A[n] = A[n] .* ( ( tenmat(X,n) * SAn' ) ./ (A[n] * SAn * SAn' ))
        end
        S = X
        for m=1:N
            S = ttm(S, A[m]', m)
        end
        S .= max.(S, eps)

        numerator = X
        denominator = S
        for m=1:N
            numerator = ttm(numerator, A[m]', m)
            denominator = ttm(denominator, A[m]'*A[m], m)
        end
        S .= S .* ( numerator ./ denominator )
    end
    return S, A
end

function update_An_LS(A, n, SAn, X)
    An = A[n] .* ( ( tenmat(X,n) * SAn' ) ./ (A[n] * SAn * SAn' ))
   return An
end

function update_An_KL(A, n, SAn, X)
    tensor_size = size(X)
    z = sum(SAn, dims=2)
    An = A[n] .* ((( tenmat(X,n) ./ ( A[n] * SAn ) ) * SAn' ) ./ ( ones(tensor_size[n]) * z' ))
    return An
end

function update_S_LS(S, A, X)
    N = ndims(S)
    numerator = X
    denominator = S
    for m=1:N
        numerator = ttm(numerator, A[m]', m)
        denominator = ttm(denominator, A[m]'*A[m], m)
    end
    S = S .* ( numerator ./ denominator )
    return S
end

function update_S_KL(S, A, X)
    N = ndims(S)
    numerator = X ./ reconst(S, A)
    denominator = ones( size(X)... )
    for m=1:N
        numerator = ttm(numerator, A[m]', m)
        denominator = ttm(denominator, A[m]', m)
    end
    S = S .* ( numerator ./ denominator )
    return S
end

function update(X, S, A, cost, max_iter=100, verbose=true, verbose_interval=20)
    N = ndims(S)

    cnt_iter = 0
    while(cnt_iter < max_iter)
        ############
        # update A #
        ############
        for n=1:N
            # get SAn
            SAn = S
            for m in [1:n-1;n+1:N]
                SAn = ttm(SAn, A[m], m)
            end
            SAn = tenmat(SAn, n)
            if cost == "LS"
                A[n] .= update_An_LS(A, n, SAn, X)
            elseif cost == "KL"
                A[n] .= update_An_KL(A, n, SAn, X)
            end
        end

        ########################
        # update Core tensor S #
        ########################
        if cost == "LS"
            S .= update_S_LS(S, A, X)
        elseif cost == "KL"
            S .= update_S_KL(S, A, X)
        end

        if verbose && (cnt_iter % verbose_interval == 0)
            Xr = reconst(S, A)
            error = cost_function(cost, X, Xr)
            println("$cnt_iter $error")
        end
        cnt_iter += 1
    end
    return S, A
end

"""
Non-Negative Tucker Decomposition
proposed by Young-Deok Kim et al.
See [original paper](https://ieeexplore.ieee.org/document/4270403)

# Aruguments
- `X` : input non-negative tensor
- `reqrank` : Target Tucker rank, array
- `init_method` : initial values, "NMF" or "random"
- `cost` : cost function, "LS" or "KL"
- `verbose` : true or false
- `pre_iter_max` : iter_max of initialization based on "NMF"
- `pre_iter_max_inloop` : iter_max of initialization of "NMF"
"""
function NTD(X, reqrank ;
        cost="LS", init_method="NMF", max_iter=500, verbose=true, verbose_interval=50,
        pre_iter_max=40, pre_iter_max_inloop=40)

    @assert ndims(X) == length(reqrank)
    tensor_size = size(X)

    #A[n] \in R^( tensor_size[n] \times reqrank[n] )
    #S \in R^(reqrank[1] \times ... \times reqrank[N])
    if init_method == "random"
        S, A = init_random(tensor_size, reqrank)
    elseif init_method == "NMF"
        S, A = init_based_NMF(X, reqrank, pre_iter_max=pre_iter_max, pre_iter_max_inloop=pre_iter_max_inloop)
    else
        error("no init method ", init_method)
    end

    S, A = update(X, S, A, cost, max_iter, verbose, verbose_interval)
    Xr = reconst(S, A)
    return S, A, Xr
end

Xrが再構成後の低ランクテンソル．mrank(Xr)でXrのランクが落ちていることを確認できる．

Julia言語 Non-Negative Tucker Decomposition の実装 [LSエラー版]

↓　こっちのほうが役に立つと思う

genkaiphd.hatenablog.com

↑　こっちのほうが役に立つと思う

研究に必要なので，非負テンソルのタッカー分解の論文をjuliaで実装した．とりあえず，コスト関数がLSエラーのものを実装した．KLダイバージェンスの場合についてはまた明日，実装する予定．

ieeexplore.ieee.org

論文中にいくつかトラップがある．

・式(17)は間違い．（まぁ，これは実装には関係ないけど．）

・Table4の S←Rand(R,m)はS←Rand(R,l)が正しいと思われる．

・Table4, Table5に無定義で出てくる行列の肩の＋は転置Tの誤植と思われる．

あと，式(22)を使うと計算効率がよくなるっぽい．

using LinearAlgebra
using TensorToolbox

function reconst(S, A)
    N = ndims(S)
    Xr = ttm(S, A[1], 1)
    for n=2:N
        Xr = ttm(Xr, A[n], n)
    end
    return Xr
end

function init_random(tensor_size, reqrank)
    A = []
    N = length(reqrank)
    for n = 1:N
        An = rand(tensor_size[n], reqrank[n])
        push!(A, An)
    end
    S = rand( reqrank... )
    return S, A
end

function init_for_NMF(X, R, eps=1.0e-10, iter_max_inloop=40)
    (m, l) = size(X)
    A = rand(m, R)
    S = rand(R, l)

    for iter = 1:iter_max_inloop
        A .= max.(X*S', eps)
        A .= A .* ((X*S') ./ (A*S*S'))
        S .= max.(A'*X, eps)
        S .= S .* ( (A'*X) ./ (A'*A*S) )
    end
    return S, A
end

function init_based_NMF(X, reqrank, eps=1.0e-10, iter_max=40)
    N = ndims(X)
    A = []
    for n=1:N
        _, An = init_for_NMF(tenmat(X, n), reqrank[n])
        push!(A, An)
    end

    S = rand( reqrank... )
    for iter = 1:iter_max
        SAns = []
        for n=1:N
            SAn = S
            for m=1:N
                if m != n
                    SAn = ttm(SAn, A[m], m)
                end
            end
            SAn = tenmat(SAn, n)
            push!(SAns, SAn)
            A[n] = A[n] .* ( ( tenmat(X,n) * SAns[n]' ) ./ (A[n] * SAns[n] * SAns[n]' ))
        end
        S = X
        for m=1:N
            S = ttm(S, A[m]', m)
        end
        S .= max.(S, eps)

        numerator = X
        denominator = S
        for m=1:N
            numerator = ttm(numerator, A[m]', m)
            denominator = ttm(denominator, A[m]'*A[m], m)
        end
        S .= S .* ( numerator ./ denominator )
    end
    return S, A
end

function update(X, S, A, max_iter=100, verbose=true, verbose_interval=20)
    N = ndims(S)

    cnt_iter = 0
    while(cnt_iter < max_iter)
        ############
        # update A #
        ############
        SAns = []
        for n=1:N
            # get SAn
            SAn = S
            for m=1:N
                if m != n
                    SAn = ttm(SAn, A[m], m)
                end
            end
            SAn = tenmat(SAn, n)
            push!(SAns, SAn)
            A[n] = A[n] .* ( ( tenmat(X,n) * SAns[n]' ) ./ (A[n] * SAns[n] * SAns[n]' ))
        end

        ########################
        # update Core tensor S #
        ########################
        numerator = X
        denominator = S
        for m=1:N
            numerator = ttm(numerator, A[m]', m)
            denominator = ttm(denominator, A[m]'*A[m], m)
        end
        S .= S .* ( numerator ./ denominator )

        if verbose && (cnt_iter % verbose_interval == 0)
            Xr = reconst(S, A)
            cost = norm(X - Xr)
            println("$cnt_iter $cost")
        end
        cnt_iter += 1
    end

    return S, A
end

"""
Non-Negative Tucker Decomposition
proposed by Young-Deok Kim et al.
See [original paper](https://ieeexplore.ieee.org/document/4270403)

# Aruguments
- `X` : input non-negative tensor
- `reqrank` : Target Tucker rank, array
- `init_method` : initial values, "NMF" or "random"
"""
function NTD(X, reqrank ; init_method="NMF", max_iter=400, verbose=true, verbose_interval=20)
    @assert ndims(X) == length(reqrank)
    tensor_size = size(X)

    #A[n] \in R^( tensor_size[n] \times reqrank[n] )
    #S \in R^(reqrank[1] \times ... \times reqrank[N])
    if init_method == "random"
        S, A = init_random(tensor_size, reqrank)
    elseif init_method == "NMF"
        S, A = init_based_NMF(X, reqrank)
    else
        error("no init method ", init_method)
    end

    S, A = update(X, S, A, max_iter, verbose, verbose_interval)
    Xr = reconst(S, A)

    return S, A, Xr
end

実装に丸一日かけてしまった．

Julia言語 行列の各要素が1.0e-10以下だったら0にする

一番手っ取り早いのはこんな感じ

A = rand(5,5)
A[A.<1.0e-10] .= 0

以下，強い人に教わったこと． 絶対値が1e-10以下なら0にする場合に真っ先に思い浮かぶのは

X[abs.(X) .≤ 1e-10] .= 0

だけど，これは遅い． やりたいこと分解して，

small2zero(x, ε=1e-10) = ifelse(abs(x) ≤ ε, zero(x), x)
A .= small2zero.(A)

とすればよい．

#Julia言語 絶対値が1e-10以下なら0にする函数を書く:

small2zero(x, ε=1e-10) = ifelse(abs(x) ≤ ε, zero(x), x)

配列Aに

A .= small2zero.(A)

または

@. A = small2zero(A)

のように適用する。

やりたいことを分割して、小さな汎用的函数の集まりで問題を解くと楽。 https://t.co/CxetvrrVOg pic.twitter.com/iB7VTtCtR3

— 黒木玄 Gen Kuroki (@genkuroki) 2021年4月8日