tucker-hosvd.jl (1379B)
1 #/usr/bin/julia 2 3 using LinearAlgebra 4 using Plots 5 using Distributions 6 using LaTeXStrings 7 using Random 8 using TSVD: tsvd # todo: implement 9 using TensorToolbox: tenmat, ttm, contract # todo:implement 10 using Plots.Measures 11 12 function my_hosvd(C, r) 13 S0 = copy(C); d = length(size(C)); n = size(C) 14 V = []; singular_values = []; errors = []; 15 for k=1:d 16 # reshape & permute S_{k-1} -> B 17 B_perm = permutedims(S0, ([[i for i=1:d if i!=k]..., k])) 18 B = reshape(B_perm, (prod([[r_i for r_i in r[1:(k-1)]]..., [i_k for i_k in n[(k+1):d]]...])..., n[k])) 19 B = convert(Matrix{Float64}, B) 20 21 # rank r_k T-SVD of B -> B^hat 22 U_hat, Sig_hat, V_hat = tsvd(B, r[k]) 23 B_hat = U_hat * Diagonal(Sig_hat) * transpose(V_hat) 24 25 # reshape & permute B^hat -> S_k 26 W_hat = B_hat * V_hat 27 W_reshape = reshape(W_hat , ([r_i for r_i in r[1:(k-1)]]..., [n_k for n_k in n[k+1:end]]..., r[k])) 28 S0 = permutedims(W_reshape , ([i for i=1:k-1]..., d, [i for i=k:d-1]...)) 29 30 31 append!(V, [V_hat]) 32 append!(singular_values, [Sig_hat]) 33 C_hat = tucker_eval(S0, V) 34 append!(errors, norm(C_hat - C)/norm(C)) 35 end 36 return V, S0, singular_values, errors 37 end 38 39 function tucker_eval(S, V) 40 d = length(V); A = copy(S) 41 for (k, V_k) in enumerate(V) 42 A = ttm(A, V_k, k) 43 end 44 return A 45 end