まぃふぇいばりっと

退職D進して機械学習やってます.Julia大好きです.勉強したことの殴り書きです.

Julia言語 L1ノルムでのテンソルランク1近似の解析解

3階のテンソルをL1ノルムの尺度で近似するランク1テンソルを見つけます,という論文を実装した.論文のタイトルにexact solutionがあるけど,L1ノルムの意味での最良1ランク近似をしてるわけではないので注意(かなり混乱した).

arxiv.org

最近,L1が流行っているのかな.著者にお願いしてuploadしてもらったpythonコードとまったく同じ結果を返したので,このjulia実装は正しいと思う.私はgreedyのバージョンしか実装していないので,もっと速くする方法も論文を読めば載ってます(実装大変そうだからあきらめたけど,pythonコードを翻訳すればいけるかも?)

github.com

using TensorToolbox
using LinearAlgebra
using Arpack

function exactL1(tensor)
    X = tenmat(tensor, 1) #mode-1 matricization
    (D, M, N) = size(tensor)
    B = vec(collect(Iterators.product(fill([-1,1],N)...)))
    numOfCandidates = 2^N

    metopt = 0
    bopt = undef
    for n = 1:numOfCandidates
        b = collect(B[n])
        Z = X * kron(b, Matrix(I, M, M))
        sigmamax = svds(Z, nsv=1)[1].S[1]
        if sigmamax > metopt
            metopt = sigmamax
            bopt = b
        end
    end

    USV = svds( X*kron(bopt, Matrix(I,M,M)), nsv=1)[1]
    uopt = USV.U[:,1]
    vopt = USV.Vt[1,:]
    return uopt, vopt, bopt, metopt
end

入力行列の復元は原論文に従って,こんな感じ.試しにランク1テンソルを入力にしてみると,エラーがとても小さくなることが確かめられた.

P = [1 2]
Q = [2 3 4]
R = [4 5 6 7]
T = dropdims( ttt(ttt(P,Q),R) )
uopt,vopt,bopt,metopt = exactL1(T)

T1 = zeros(2,3,4)
for i = 1 : 4
    T1[:,:,i] = ( uopt * uopt' ) * T[:, :, i] * ( vopt * vopt' )
end

println( norm(T-T1) )