main.jl (10779B)
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 include("mps-tt.jl") 13 include("tucker-hosvd.jl") 14 include("unf-aprx.jl") 15 include("functions.jl") 16 17 function main() 18 # Exercise 2 19 d = 4; 20 n = [20+k for k=1:d]; 21 r = [2*k for k=1:d]; 22 V = [rand(Uniform(-1, 1), n[k], r[k]) for k=1:d]; 23 S = rand(Uniform(-1, 1), r...); 24 C = tucker_eval(S, V); 25 Vs, S0, vals, errors = my_hosvd(C, r); 26 27 p = plot(errors, 28 ylabel="error", 29 label=L"$\frac{\|\hat{A}_k - A\|_F}{\|A\|_F}$", 30 xlabel=L"k", 31 lw = 1, 32 xticks=collect(1:d), 33 markershape=:circle, 34 legend=:topleft, 35 margin=5mm, 36 dpi=300, 37 title="HOSVD of Uniform Tensor") 38 savefig(p, "./plots/hosvd-uniform-error.png") 39 40 # Exercise 3 41 d = 4; n = 51; 42 A = [f([t1, t2, t3, t4]) for t1 in t(n), t2 in t(n), t3 in t(n), t4 in t(n)]; 43 B = [g([t1, t2, t3, t4]) for t1 in t(n), t2 in t(n), t3 in t(n), t4 in t(n)]; 44 45 # Exercise 3a 46 σ_unf_a = σ_unfold(A, d); 47 σ_unf_b = σ_unfold(B, d); 48 p1 = plot(margin=5mm) 49 p2 = plot(margin=5mm) 50 for k=1:d 51 plot!(p1, 52 σ_unf_a[k], 53 ylabel=L"$\log(\Sigma_{A_k})$", 54 label=L"$\Sigma_{A_%$k}$", 55 xlabel=L"n", 56 lw = 1, 57 yaxis=:log, 58 markershape=:circle, 59 dpi=300, 60 title="Unfolding of A") 61 plot!(p2, 62 σ_unf_b[k], 63 ylabel=L"$\log(\Sigma_{B_k})$", 64 label=L"$\Sigma_{_B%$k}$", 65 xlabel=L"n", 66 lw = 1, 67 dpi=300, 68 markershape=:circle, 69 yaxis=:log, 70 title="Unfolding of B") 71 end 72 savefig(p1, "./plots/singular-dec-A.png") 73 savefig(p2, "./plots/singular-dec-B.png") 74 75 _s = [1/(10^(j*2)) for j=1:5] # computer not goode enough for j = 6 76 77 r_jk_a, ϵ_jk_a, σ_jk_a = rank_approx(A, tucker_unfold); 78 r_jk_b, ϵ_jk_b, σ_jk_b = rank_approx(B, tucker_unfold); 79 p1 = plot(size=[700, 350], margin=5mm, dpi=300) 80 p2 = plot(size=[700, 350], margin=5mm, dpi=300) 81 for k=1:d 82 plot!(p1, 83 1 ./ ϵ_s, 84 r_jk_a[:, k], 85 ylabel=L"$r_{jk}$", 86 label=L"A-$\varepsilon_{j%$k}$", 87 xlabel=L"\varepsilon_{j}^{-1}", 88 lw = 2, 89 markershape=:circle, 90 xaxis=:log10, 91 legend=:topleft, 92 title=L"Tensor $f(x_1, x_2, x_3, x_4)$") 93 plot!(p2, 94 1 ./ ϵ_s, 95 r_jk_b[:, k], 96 ylabel=L"$r_{jk}$", 97 label=L"B-$\varepsilon_{j%$k}$", 98 xlabel=L"\varepsilon_{j}^{-1}", 99 lw = 2, 100 markershape=:circle, 101 xaxis=:log10, 102 legend=:topleft, 103 title=L"Tensor $g(x_1, x_2, x_3, x_4)$") 104 end 105 savefig(p1, "./plots/hosvd-error-A.png") 106 savefig(p2, "./plots/hosvd-error-B.png") 107 108 # Exercise 3b 109 Nj_hosvd_a = []; ϵ_hosvd_a = []; σ_hosvd_a = []; 110 Nj_hosvd_b = []; ϵ_hosvd_b = []; σ_hosvd_b = []; 111 112 println("Checking validity of error analysis for HOSVD") 113 for j=1:length(ϵ_s) 114 S0_a, Vs_a, vals_a, errors_a = my_hosvd(A, r_jk_a[j, :]) 115 N_vals = 0 116 for V in Vs_a 117 N_vals += length(V) 118 end 119 append!(Nj_hosvd_a, N_vals + length(S0_a)) 120 append!(ϵ_hosvd_a, errors_a[end]) 121 append!(σ_hosvd_a, [vals_a]) 122 println("ϵ_hosvd_a[j] <= |ϵ_jk_a|_f : $(ϵ_hosvd_a[j] <= norm(ϵ_jk_a)) ϵ_hosvd_a[j] = $(ϵ_hosvd_a[j])") 123 124 S0_b, Vs_b, vals_b, errors_b = my_hosvd(B, r_jk_b[j, :]) 125 N_vals = 0 126 for V in Vs_b 127 N_vals += length(V) 128 end 129 append!(Nj_hosvd_b, N_vals + length(S0_b)) 130 append!(ϵ_hosvd_b, errors_b[end]) 131 append!(σ_hosvd_b, [vals_b]) 132 println("ϵ_hosvd_b[j] <= |ϵ_jk_b|_f : $(ϵ_hosvd_b[j] <= norm(ϵ_jk_b)) ϵ_hosvd_b[j] = $(ϵ_hosvd_b[j])") 133 end 134 135 # Exercise 3b 136 j = 5 137 p1 = plot(dpi=300, size=[800, 350], margin=5mm) 138 p2 = plot(dpi=300, size=[800, 350], margin=5mm) 139 for k=1:d 140 plot!(p1, σ_unf_a[k][1:r_jk_a[j, k]] ./ σ_hosvd_a[j][k], 141 ylabel=L"\frac{\sigma^{A_k}_\alpha}{\sigma^{hosvd}_{k\alpha}}", 142 label=L"k=%$k", 143 xlabel=L"j", 144 xticks=collect(1:length(r_jk_a)), 145 yticks=[1, 1.001], 146 ylim=[1, 1.001], 147 lw = 2, 148 markershape=:circle, 149 legend=:topleft, 150 title=L"$f(x_1, x_2, x_3, x_4)$") 151 plot!(p2, σ_unf_b[k][1:r_jk_b[j, k]] ./ σ_hosvd_b[j][k], 152 ylabel=L"\frac{\sigma^{B_k}_\alpha}{\sigma^{hosvd}_{k\alpha}}", 153 label=L"k=%$k", 154 xlabel=L"j", 155 yticks=[1, 1.001], 156 ylim=[1, 1.001], 157 lw = 2, 158 markershape=:circle, 159 legend=:topleft, 160 title=L"$g(x_1, x_2, x_3, x_4)$") 161 end 162 savefig(p1, "./plots/hosvd-sigmaratio-a.png") 163 savefig(p2, "./plots/hosvd-sigmaratio-b.png") 164 165 # Exercise 3d 166 167 p1 = plot(size=[700, 350], margin=5mm, dpi=300) 168 p2 = plot(size=[700, 350], margin=5mm, dpi=300) 169 plot!(p1, 170 1 ./ ϵ_s, 171 Nj_hosvd_a, 172 ylabel=L"$N_{j}$", 173 label=L"A-$N_{j}$", 174 xlabel=L"\varepsilon_{j}^{-1}", 175 lw = 2, 176 markershape=:circle, 177 xaxis=:log10, 178 legend=:topleft, 179 title=L"Tucker of $f(x_1, x_2, x_3, x_4)$") 180 plot!(p2, 181 1 ./ ϵ_s, 182 Nj_hosvd_b, 183 ylabel=L"$N_{j}$", 184 label=L"B-$N_{j}$", 185 xlabel=L"\varepsilon_{j}^{-1}", 186 lw = 2, 187 markershape=:circle, 188 xaxis=:log10, 189 legend=:topleft, 190 title=L"Tucker of $g(x_1, x_2, x_3, x_4)$") 191 savefig(p1, "./plots/hosvd-Nj-a.png") 192 savefig(p2, "./plots/hosvd-Nj-b.png") 193 194 # exercise 6 195 d = 4; 196 n = [20+k for k=1:d]; 197 r = [[2*k for k=1:(d-1)]..., 1]; 198 r_new = [1, r...] 199 V = [rand(Uniform(-1, 1), (r_new[i], n[i], r_new[i+1])) for i=1:d]; 200 D = ttmps_eval(V, n) 201 C, singular_val, errors= tt_svd(D, n, r, d); 202 p = plot(errors, 203 ylabel="error", 204 label=L"$\frac{\|\hat{A}_k - A\|_f}{\|A\|_f}$", 205 xlabel=L"k", 206 lw = 1, 207 xticks=collect(1:d), 208 markershape=:circle, 209 legend=:topleft, 210 margin=5mm, 211 dpi=300, 212 title="TT-SVD of uniform tensor") 213 savefig(p, "./plots/ttsvd-uniform-error.png") 214 215 # exercise 7 216 217 d = 4; n = 51; 218 ndim = [n for i=1:d] 219 A = [f([t1, t2, t3, t4]) for t1 in t(n), t2 in t(n), t3 in t(n), t4 in t(n)]; 220 B = [g([t1, t2, t3, t4]) for t1 in t(n), t2 in t(n), t3 in t(n), t4 in t(n)]; 221 222 ## Exercise 7a 223 224 r_jk_a, ϵ_jk_a, σ_jk_a = rank_approx(A, ttmps_unfold); 225 r_jk_b, ϵ_jk_b, σ_jk_b = rank_approx(B, ttmps_unfold); 226 p1 = plot(size=[700, 350], margin=5mm, dpi=300) 227 p2 = plot(size=[700, 350], margin=5mm, dpi=300) 228 for k=1:d 229 plot!(p1, 230 1 ./ ϵ_s, 231 r_jk_a[:, k], 232 ylabel=L"$r_{jk}$", 233 label=L"A-$\varepsilon_{j%$k}$", 234 xlabel=L"\varepsilon_{j}^{-1}", 235 lw = 2, 236 markershape=:circle, 237 xaxis=:log10, 238 legend=:topleft, 239 title=L"Tensor $f(x_1, x_2, x_3, x_4)$") 240 plot!(p2, 241 1 ./ ϵ_s, 242 r_jk_b[:, k], 243 ylabel=L"$r_{jk}$", 244 label=L"B-$\varepsilon_{j%$k}$", 245 xlabel=L"\varepsilon_{j}^{-1}", 246 lw = 2, 247 markershape=:circle, 248 xaxis=:log10, 249 legend=:topleft, 250 title=L"Tensor $g(x_1, x_2, x_3, x_4)$") 251 end 252 savefig(p1, "./plots/tucker-error-A.png") 253 savefig(p2, "./plots/tucker-error-B.png") 254 255 # exercise 7b 256 Nj_ttsvd_a = []; ϵ_ttsvd_a = []; σ_ttsvd_a = []; 257 Nj_ttsvd_b = []; ϵ_ttsvd_b = []; σ_ttsvd_b = []; 258 259 println("Checking validity of error analysis for TT-SVD") 260 for j=1:size(ϵ_jk_a, 1) 261 Vs_a, vals_a, errors_a = tt_svd(A, ndim, r_jk_a[j, :], d); 262 N_vals = 0 263 for V in Vs_a 264 N_vals += length(V) 265 end 266 append!(Nj_ttsvd_a, N_vals) 267 append!(ϵ_ttsvd_a, errors_a[end]) 268 append!(σ_ttsvd_a, [vals_a]) 269 println("ϵ_ttsvd_a[j] <= |ϵ_jk_a|_f : $(ϵ_ttsvd_a[j] <= norm(ϵ_jk_a)) ϵ_ttsvd_a[j] = $(ϵ_ttsvd_a[j])") 270 271 Vs_b, vals_b, errors_b = tt_svd(B, ndim, r_jk_b[j, :], d); 272 N_vals = 0 273 for V in Vs_b 274 N_vals += length(V) 275 end 276 append!(Nj_ttsvd_b, N_vals) 277 append!(ϵ_ttsvd_b, errors_b[end]) 278 append!(σ_ttsvd_b, [vals_b]) 279 println("ϵ_ttsvd_b[j] <= |ϵ_jk_b|_f : $(ϵ_ttsvd_b[j] <= norm(ϵ_jk_b)) ϵ_ttsvd_b[j] = $(ϵ_ttsvd_b[j])") 280 end 281 282 j = 5 283 p1 = plot(dpi=300, size=[800, 350], margin=5mm) 284 p2 = plot(dpi=300, size=[800, 350], margin=5mm) 285 for k=1:d-1 286 plot!(p1, σ_unf_a[k][1:r_jk_a[j, k]] ./ σ_ttsvd_a[j][k], 287 ylabel=L"\frac{\sigma^{A_k}_\alpha}{\sigma^{ttsvd}_{k\alpha}}", 288 label=L"k=%$k", 289 xlabel=L"j", 290 xticks=collect(1:length(r_jk_a)), 291 yticks=[1, 1000], 292 # ylim=[1, 1000], 293 lw = 2, 294 markershape=:circle, 295 legend=:topleft, 296 title=L"$f(x_1, x_2, x_3, x_4)$") 297 plot!(p2, σ_unf_b[k][1:r_jk_b[j, k]] ./ σ_ttsvd_b[j][k], 298 ylabel=L"\frac{\sigma^{B_k}_\alpha}{\sigma^{ttsvd}_{k\alpha}}", 299 label=L"k=%$k", 300 xlabel=L"j", 301 yticks=[1, 1000], 302 # ylim=[1, 1000], 303 lw = 2, 304 markershape=:circle, 305 legend=:topleft, 306 title=L"$g(x_1, x_2, x_3, x_4)$") 307 end 308 savefig(p1, "./plots/ttsvd-sigmaratio-a.png") 309 savefig(p2, "./plots/ttsvd-sigmaratio-b.png") 310 311 312 # Exercise 7d 313 p1 = plot(size=[700, 350], margin=5mm, dpi=300) 314 p2 = plot(size=[700, 350], margin=5mm, dpi=300) 315 plot!(p1, 316 1 ./ ϵ_s, 317 Nj_ttsvd_a, 318 ylabel=L"$N_{j}$", 319 label=L"A-$N_{j}$", 320 xlabel=L"\varepsilon_{j}^{-1}", 321 lw = 2, 322 markershape=:circle, 323 xaxis=:log10, 324 legend=:topleft, 325 title=L"MPS-TT of $f(x_1, x_2, x_3, x_4)$") 326 plot!(p2, 327 1 ./ ϵ_s, 328 Nj_ttsvd_b, 329 ylabel=L"$N_{j}$", 330 label=L"B-$N_{j}$", 331 xlabel=L"\varepsilon_{j}^{-1}", 332 lw = 2, 333 markershape=:circle, 334 xaxis=:log10, 335 legend=:topleft, 336 title=L"MPS-TT of $g(x_1, x_2, x_3, x_4)$") 337 savefig(p1, "./plots/ttsvd-Nj-a.png") 338 savefig(p2, "./plots/ttsvd-Nj-b.png") 339 340 end 341 342 main()