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

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

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

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

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) )