まぃふぇいばりっと

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

Julia言語 高速3階テンソル1ランク近似手法SeROAPの実装

3階のテンソルTをST-HOSVD, T-HOSVDより正確に高速に1ランク近似できる手法をJulia言語で実装した. あ,ここでいう階っていうのは,rank(T)じゃなくて,ndims(T)ね.3-waysとか,3-orderってこと.日本では,ndimsを階っていうので,(CP)ランクと混ざって良く分からなくなる.この手法は3階テンソル以外にも使えるけど,ST-HOSVDやT-HOSVDより,より近似をする理論的保証があるのは3階のときのみ. ちょっとめんどくさかったので,入力行列Tの型をI1×I2×...×Inとしたとき,I1 >= I2 >= … を仮定している.

アルゴリズム詳細はこの論文のアルゴリズム2を参照.

arxiv.org

実装

using TensorToolbox
using LinearAlgebra
using Arpack
using TensorDecompositions

function SeROAP(T)
    # arxiv.org/pdf/1508.05273.pdf
    # input_tensor_shape = [I1, I2, ... , In]
    #  I1 >= I2 >= ... >= In have to be satisfied ?
    # input is any complex tensor
    # output is rank-1 complex tensor
    # cost function is L2
    # in case of 3-ways tensor, SeROAP is better than HOSVD

    N = ndims(T)
    input_tensor_shape = size(T)

    Vs = []
    vs = []
    V = tenmat(T,1)
    V0 = copy(V)
    for n=1:N-2
        push!(vs, svds( V ; nsv=1 )[1].V[:,1])
        push!(Vs, reshape(vs[n], input_tensor_shape[n+1], prod(input_tensor_shape[n+2:N])))
        V = Vs[n]
    end
    u = svds( V ; nsv=1 )[1].U[:, 1]
    v = svds( V ; nsv=1 )[1].V[:, 1]
    w = kron( conj(v), u)

    X_n = undef
    for n = N-2 : -1 : 1
        if n != 1
            X_n = ( Vs[n-1] * w ) * w'
        else
            X_n = ( V0 * w ) * w'
        end
        w = vec(X_n)
    end

    X = matten(X_n, 1, collect(input_tensor_shape) )
    return X
end

研究には使わないけど,一緒にDCPDを実装しておこうかな.