SVDによるテンソルリング分解を実装しました．テンソルリング分解についてはこちらの論文を参照．テンソルリング分解では，ランクはベクトルです．このベクトルの各要素には自然数がはいっていて，要素数は入力テンソルのオーダーに一致します．
こちらにも実装の詳細があって，結構助かりました．
オリジナルの論文では，リングランクを指定するのではなくて，許容する誤差εを入力に与えています．自分の研究ではリングランクを決め打ちする実装も必要になりそうでした．そこで，TR()の二つ目の引数にベクトルを渡したら，ランクを決め打ちで分解，floatを渡したら誤差ε以下で分解，という実装にしました．多重ディスパッチ最高．
自然数の約数を求めるアルゴリズムと，誤差を指定してtruncated SVD が必要だったので，こちらを先に用意しておきます．約数を全部求めるアルゴリズムはお気楽 Julia プログラミング超入門から拝借しました．ありがとうございます。
""" Provided by http://www.nct9.ne.jp/m_hiroi/light/juliaa08.htmll """ function divisor(n) # pq の約数を求める divisor_sub(p, q) = [p ^ i for i = 0 : q] xs = factorization(n) ys = divisor_sub(xs[1]...) for p = xs[2 : end] ys = [x * y for x = divisor_sub(p...) for y = ys] end sort(ys) end function factorSub(n, m) c = zero(n) while n % m == 0 c += 1 n = div(n, m) end c, n end function factorization(n) x::typeof(n) = 2 xs = Pair{typeof(n), typeof(n)}[] c, m = factorSub(n, x) if c > 0 push!(xs, x => c) end x = 3 while x * x <= m c, m = factorSub(m, x) if c > 0 push!(xs, x => c) end x += 2 end if m > 1 push!(xs, m => one(n)) end xs end
truncated SVD については，いままでいろいろ実装してきましたが，全てランクを先に決め打ちするものでした．今回はランクを決め打ちするものと，許容誤差をあたえておく低ランク近似の両方の関数を用意しておきます．
function lowrank_app_th(M,δ) m = rank(M) U,Σ,V = svd(M) sorted_Σ = sort(Σ) err = cumsum(sorted_Σ[begin:end-1].^2) reverse!(err) r = m for l = 1:m-1 if err[l] < δ r = l break end end Ur = U[:,1:r] Σr = Σ[1:r] Vr = V[:,1:r] # reconstraction is Ur*diagm(Σr)*Vr' return Ur,diagm(Σr),Vr,r end function lowrank_app(M,r::Int) U,Σ,V = svd(M) Ur = U[:,1:r] Σr = Σ[1:r] Vr = V[:,1:r] # reconstraction is Ur*diagm(Σr)*Vr' return Ur,diagm(Σr),Vr end
ここからが本体．分解したコアテンソルから元のテンソルを再構成する関数reconst()も用意しておきます．
どうもreconst()では結構時間がかかるので，もっとよい実装があるはず...。
using LinearAlgebra using Glob using Plots using Images using Statistics """ Tensor Ring Decomposition based on SVD proposed by [Q Zhao (2016)](https://arxiv.org/abs/1606.05535) [Matlab code](https://github.com/oscarmickelin/tensor-ring-decomposition) is also available, # Aruguments - 'T' : input tensor - 'r' : target rank, which should be vector example: T = randn(20,20,30) r = (5,3,10) ring_cores = TR(T,r) size.(ring_cores) # 3-element Vector{Tuple{Int64, Int64, Int64}}: # (5, 20, 3) # (3, 20, 10) # (10, 30, 5) # We can obtain low-ring rank tensor by 'reconst' Tr = reconst(ring_cores) """ function TR(T,r::Vector) n = size(T) d = length(r) @assert r[1]*r[2] <= n[1] "r1*r2 should be smaller than n[1]" @assert d == ndims(T) "the length of ranks should be same as the dims of tensor" r0 = r[1]*r[2] T1 = reshape(T,(n[1],prod(n[2:end]))) U,Σ,V = lowrank_app(T1,r0) G = [] G1 = permutedims(reshape(U,(n[1],r[1],r[2])),[2,1,3]) C = permutedims(reshape(Σ*V',(r[1],r[2],prod(n[2:d]))),[2,3,1]) push!(G,G1) for k=2:d-1 C = reshape(C,(r[k]*n[k],prod(n[k+1:end])*r[1])) U,Σ,V = lowrank_app(C,r[k+1]) Gk = reshape(U,(r[k],n[k],r[k+1])) push!(G,Gk) C = reshape(Σ*V',(r[k+1],:,r[1])) end Gd = reshape(C,(r[d],n[d],r[1])) push!(G,Gd) return G end """ Tensor ring decomposition with prescribed error value ε. The reconstraction error becomes smaller than ε. See more details in the [paper](https://arxiv.org/abs/1606.05535) # Aruguments - 'T' : input tensor - 'ε' : prescribed error, which should be float64 example: Tr = reconst(ring_cores) T = rand(10,20,7) ε = 0.3 G,r = TR(T,ε); Tre = reconst(G) @show norm(T - Tre)/norm(T) < ε # true Gnew = TR(T,r); for l = 1:length(r) @show Gnew[l] ≈ G[l] # true end """ function TR(T,ε::Float64) n = size(T) d = ndims(T) r = zeros(Int,d) G = [] δ1 = sqrt(2/d)*ε*norm(T) δ = δ1/sqrt(2) T1 = reshape(T,(n[1],prod(n[2:end]))) U,Σ,V, r0 = lowrank_app_th(T1,δ1) factors = divisor(r0) r[1] = factors[floor(Int,length(factors)/2)] r[2] = Int( r0 / r[1] ) G1 = permutedims(reshape(U,(n[1],r[1],r[2])),[2,1,3]) C = permutedims(reshape(Σ*V',(r[1],r[2],prod(n[2:d]))),[2,3,1]) push!(G,G1) for k=2:d-1 C = reshape(C,(r[k]*n[k],prod(n[k+1:end])*r[1])) U,Σ,V, rnew = lowrank_app_th(C,δ) r[k+1] = rnew Gk = reshape(U,(r[k],n[k],r[k+1])) push!(G,Gk) C = reshape(Σ*V',(r[k+1],:,r[1])) end Gd = reshape(C,(:,n[d],r[1])) r[d] = size(Gd)[1] push!(G,Gd) return G, r end """ Reconstract the original tensor from ring cores. # Arguments - 'G' : cores tensors """ function reconst(G) d = length(G) I = zeros(Int,d) for n = 1:d _, In, _ = size(G[n]) I[n] = In end reconst_T = zeros(I...) for idx in CartesianIndices((tuple(I...))) reconst_T[idx] = tr(prod([G[n][:,idx[n],:] for n=1:d])) end return reconst_T end
テスト
適当なテンソルを用意して，テンソルリングを大きくしていくとちゃんと誤差が小さくなることを確かめてみます．
T = rand(25,25,25) for r in [1,2,3,4,5] G = TR(T,[r,r,r]) Tre = reconst(G) @show norm(T-Tre) end # norm(T - Tre) = 36.07195656407602 # norm(T - Tre) = 35.709229547710216 # norm(T - Tre) = 35.17180322846426 # norm(T - Tre) = 34.43064761988541 # norm(T - Tre) = 33.56737731452813
問題なさそうです．
決め打ちした誤差以下になっているか確かめてみます．
for ε in [10.3, 9.4, 3.3, 1.2, 0.8, 0.01, 0.001] G, r = TR(T,ε); Tre = reconst(G); err = norm(T - Tre)/norm(Tre); @show err < ε end
全部trueになります．問題なさそうです．
最後に適当に画像を再構成してみます．ORL顔画像データセットをダウンロードして，/dataに置いておきます．410枚のグレースケール顔画像を32×27にリサイズして，まとめた32×27×410テンソルTをつくります．とりあえず１枚目の顔画像T[:,:,1]をプロットさせてみます．
names = glob("data/*") L = length(names) #410 H = 32 W = 27 T = zeros(H,W,L) l = 1 for name in names img = load(name) img = Gray.(img) img = imresize(img, (H, W)) img_mtx = convert(Array{Float64}, img) T[:,:,l] = img_mtx l += 1 end plot(Gray.(T[:,:,1]))
このテンソルTを分解して，ε=0.2で再構成して，1枚目の顔画像を復元した様子がこちら．
G, r =TR(T,0.2); T_reconst = reconst(G); plot(Gray.(T_reconst[:,:,1])) @show r #(4,7,162)
このときランクは(4,7,162)でした．
ランクを小さくすると当然再構成が雑になっていきます．順に第三リングランクを50と10にしてみた様子です．
G=TR(T,[4,7,50]); T_reconst = reconst(G); plot(Gray.(T_reconst[:,:,1])) G=TR(T,[4,7,50]); T_reconst = reconst(G); plot(Gray.(T_reconst[:,:,1]))
reconst関数のもっと速い実装
Lirimyさん(@LirimyDh)に次のような高速な実装を教わりました．すごい速くなりました．ありがとうございます．以下，コードの一部をこのページから引用します．私もこういうコードがかけるようになりたいものです．
# ３次元配列を行列の配列に変換する "Just specify that [an element of] the argument is an array of matrices" struct MatrixArray end """ R = matrixarray(S::AbstractArray{T,3}) `R[k][i, j] = S[i, k, j]` """ matrixarray(S) = [SMatrix{size(S[:, i, :])...}(S[:, i, :]) for i in axes(S, 2)] #matrixarray(S) = [S[:, i, :] for i in axes(S, 2)] # without StaticArrays # 行列積の最終段を省き、トレースを直接求める reconst(G) = reconst(MatrixArray(), matrixarray.(G)) reconst(::MatrixArray, G) = _reconst(G...) # Gs を真ん中で分けないと遅くなる function _reconst(Gs...) h = length(Gs) ÷ 2 _reconst(reduce(eachprod, Gs[begin:h]), reduce(eachprod, Gs[h+1:end])) end #_reconst(G1, G2, Gs...) = _reconst(eachprod(G1, G2), Gs...) # 遅い _reconst(G1, G2) = trprod.(G1, expanddim(G2, G1)) """ trprod(A, B) Returns `tr(A * B)` """ trprod(A, B) = dot(vec(A'), vec(B)) """ C = eachprod(A, B) `C[i, j, k] = A[i, j] * B[k]` `A, B, C :: Array{Matrix}` """ eachprod(A, B) = A .* expanddim(B, A) """ Bx = expanddim(B, A) `Bx = reshape(B, (1, 1, 1, m, n))` where `ndims(A) == 3`, `size(B) == (m, n)` """ expanddim(B, A) = reshape(B, (ntuple(_ -> 1, ndims(A))..., size(B)...))
