commit 26bcda26277a4ca3877e9915a319663c9a80a49e
parent 4dcf6389aefd295070c3a96934287983f6929387
Author: miksa234 <milutin@popovic.xyz>
Date: Wed, 15 Dec 2021 10:05:20 +0100
code as05 done
Diffstat:
17 files changed, 584 insertions(+), 171 deletions(-)
diff --git a/assingments/hw05_published_2021-12-09.pdf b/assingments/hw05_published_2021-12-09.pdf
Binary files differ.
diff --git a/sesh4/src/main.jl b/sesh4/src/main.jl
@@ -16,66 +16,213 @@ include("functions.jl")
function main()
# Exercise 2
- d = 4;
- n = [20+k for k=1:d];
- r = [2*k for k=1:d];
- V = [rand(Uniform(-1, 1), n[k], r[k]) for k=1:d];
- S = rand(Uniform(-1, 1), r...);
- C = tucker_eval(S, V);
- Vs, S0, vals, errors = my_hosvd(C, r);
+ #d = 4;
+ #n = [20+k for k=1:d];
+ #r = [2*k for k=1:d];
+ #V = [rand(Uniform(-1, 1), n[k], r[k]) for k=1:d];
+ #S = rand(Uniform(-1, 1), r...);
+ #C = tucker_eval(S, V);
+ #Vs, S0, vals, errors = my_hosvd(C, r);
- p = plot(errors,
- ylabel="error",
- label=L"$\frac{\|\hat{A}_k - A\|_F}{\|A\|_F}$",
- xlabel=L"k",
- lw = 1,
- xticks=collect(1:d),
- markershape=:circle,
- legend=:topleft,
- margin=5mm,
- dpi=300,
- title="HOSVD of Uniform Tensor")
- savefig(p, "./plots/hosvd-uniform-error.png")
+ #p = plot(errors,
+ # ylabel="error",
+ # label=L"$\frac{\|\hat{A}_k - A\|_F}{\|A\|_F}$",
+ # xlabel=L"k",
+ # lw = 1,
+ # xticks=collect(1:d),
+ # markershape=:circle,
+ # legend=:topleft,
+ # margin=5mm,
+ # dpi=300,
+ # title="HOSVD of Uniform Tensor")
+ #savefig(p, "./plots/hosvd-uniform-error.png")
+
+ ## Exercise 3
+ #d = 4; n = 51;
+ #A = [f([t1, t2, t3, t4]) for t1 in t(n), t2 in t(n), t3 in t(n), t4 in t(n)];
+ #B = [g([t1, t2, t3, t4]) for t1 in t(n), t2 in t(n), t3 in t(n), t4 in t(n)];
+
+ ## Exercise 3a
+ #σ_unf_a = σ_unfold(A, d);
+ #σ_unf_b = σ_unfold(B, d);
+ #p1 = plot(margin=5mm)
+ #p2 = plot(margin=5mm)
+ #for k=1:d
+ # plot!(p1,
+ # σ_unf_a[k],
+ # ylabel=L"$\log(\Sigma_{A_k})$",
+ # label=L"$\Sigma_{A_%$k}$",
+ # xlabel=L"n",
+ # lw = 1,
+ # yaxis=:log,
+ # markershape=:circle,
+ # dpi=300,
+ # title="Unfolding of A")
+ # plot!(p2,
+ # σ_unf_b[k],
+ # ylabel=L"$\log(\Sigma_{B_k})$",
+ # label=L"$\Sigma_{_B%$k}$",
+ # xlabel=L"n",
+ # lw = 1,
+ # dpi=300,
+ # markershape=:circle,
+ # yaxis=:log,
+ # title="Unfolding of B")
+ #end
+ #savefig(p1, "./plots/singular-dec-A.png")
+ #savefig(p2, "./plots/singular-dec-B.png")
+
+ ϵ_s = [1/(10^(j*2)) for j=1:5] # computer not goode enough for j = 6
+
+ #r_jk_a, ϵ_jk_a, σ_jk_a = rank_approx(A, tucker_unfold);
+ #r_jk_b, ϵ_jk_b, σ_jk_b = rank_approx(B, tucker_unfold);
+ #p1 = plot(size=[700, 350], margin=5mm, dpi=300)
+ #p2 = plot(size=[700, 350], margin=5mm, dpi=300)
+ #for k=1:d
+ # plot!(p1,
+ # 1 ./ ϵ_s,
+ # r_jk_a[:, k],
+ # ylabel=L"$r_{jk}$",
+ # label=L"A-$\varepsilon_{j%$k}$",
+ # xlabel=L"\varepsilon_{j}^{-1}",
+ # lw = 2,
+ # markershape=:circle,
+ # xaxis=:log10,
+ # legend=:topleft,
+ # title=L"Tensor $f(x_1, x_2, x_3, x_4)$")
+ # plot!(p2,
+ # 1 ./ ϵ_s,
+ # r_jk_b[:, k],
+ # ylabel=L"$r_{jk}$",
+ # label=L"B-$\varepsilon_{j%$k}$",
+ # xlabel=L"\varepsilon_{j}^{-1}",
+ # lw = 2,
+ # markershape=:circle,
+ # xaxis=:log10,
+ # legend=:topleft,
+ # title=L"Tensor $g(x_1, x_2, x_3, x_4)$")
+ #end
+ #savefig(p1, "./plots/hosvd-error-A.png")
+ #savefig(p2, "./plots/hosvd-error-B.png")
+
+ ## Exercise 3b
+ #Nj_hosvd_a = []; ϵ_hosvd_a = []; σ_hosvd_a = [];
+ #Nj_hosvd_b = []; ϵ_hosvd_b = []; σ_hosvd_b = [];
+
+ #println("Checking validity of error analysis for HOSVD")
+ #for j=1:length(ϵ_s)
+ # S0_a, Vs_a, vals_a, errors_a = my_hosvd(A, r_jk_a[j, :])
+ # N_vals = 0
+ # for V in Vs_a
+ # N_vals += length(V)
+ # end
+ # append!(Nj_hosvd_a, N_vals + length(S0_a))
+ # append!(ϵ_hosvd_a, errors_a[end])
+ # append!(σ_hosvd_a, [vals_a])
+ # println("ϵ_hosvd_a[j] <= |ϵ_jk_a|_f : $(ϵ_hosvd_a[j] <= norm(ϵ_jk_a)) ϵ_hosvd_a[j] = $(ϵ_hosvd_a[j])")
+
+ # S0_b, Vs_b, vals_b, errors_b = my_hosvd(B, r_jk_b[j, :])
+ # N_vals = 0
+ # for V in Vs_b
+ # N_vals += length(V)
+ # end
+ # append!(Nj_hosvd_b, N_vals + length(S0_b))
+ # append!(ϵ_hosvd_b, errors_b[end])
+ # append!(σ_hosvd_b, [vals_b])
+ # println("ϵ_hosvd_b[j] <= |ϵ_jk_b|_f : $(ϵ_hosvd_b[j] <= norm(ϵ_jk_b)) ϵ_hosvd_b[j] = $(ϵ_hosvd_b[j])")
+ #end
+
+ ## Exercise 3b
+ #j = 5
+ #p1 = plot(dpi=300, size=[800, 350], margin=5mm)
+ #p2 = plot(dpi=300, size=[800, 350], margin=5mm)
+ #for k=1:d
+ # plot!(p1, σ_unf_a[k][1:r_jk_a[j, k]] ./ σ_hosvd_a[j][k],
+ # ylabel=L"\frac{\sigma^{A_k}_\alpha}{\sigma^{hosvd}_{k\alpha}}",
+ # label=L"k=%$k",
+ # xlabel=L"j",
+ # xticks=collect(1:length(r_jk_a)),
+ # yticks=[1, 1.001],
+ # ylim=[1, 1.001],
+ # lw = 2,
+ # markershape=:circle,
+ # legend=:topleft,
+ # title=L"$f(x_1, x_2, x_3, x_4)$")
+ # plot!(p2, σ_unf_b[k][1:r_jk_b[j, k]] ./ σ_hosvd_b[j][k],
+ # ylabel=L"\frac{\sigma^{B_k}_\alpha}{\sigma^{hosvd}_{k\alpha}}",
+ # label=L"k=%$k",
+ # xlabel=L"j",
+ # yticks=[1, 1.001],
+ # ylim=[1, 1.001],
+ # lw = 2,
+ # markershape=:circle,
+ # legend=:topleft,
+ # title=L"$g(x_1, x_2, x_3, x_4)$")
+ #end
+ #savefig(p1, "./plots/hosvd-sigmaratio-a.png")
+ #savefig(p2, "./plots/hosvd-sigmaratio-b.png")
+
+ ## Exercise 3d
+
+ #p1 = plot(size=[700, 350], margin=5mm, dpi=300)
+ #p2 = plot(size=[700, 350], margin=5mm, dpi=300)
+ #plot!(p1,
+ # 1 ./ ϵ_s,
+ # Nj_hosvd_a,
+ # ylabel=L"$N_{j}$",
+ # label=L"A-$N_{j}$",
+ # xlabel=L"\varepsilon_{j}^{-1}",
+ # lw = 2,
+ # markershape=:circle,
+ # xaxis=:log10,
+ # legend=:topleft,
+ # title=L"Tucker of $f(x_1, x_2, x_3, x_4)$")
+ #plot!(p2,
+ # 1 ./ ϵ_s,
+ # Nj_hosvd_b,
+ # ylabel=L"$N_{j}$",
+ # label=L"B-$N_{j}$",
+ # xlabel=L"\varepsilon_{j}^{-1}",
+ # lw = 2,
+ # markershape=:circle,
+ # xaxis=:log10,
+ # legend=:topleft,
+ # title=L"Tucker of $g(x_1, x_2, x_3, x_4)$")
+ #savefig(p1, "./plots/hosvd-Nj-a.png")
+ #savefig(p2, "./plots/hosvd-Nj-b.png")
+
+ ## exercise 6
+ #d = 4;
+ #n = [20+k for k=1:d];
+ #r = [[2*k for k=1:(d-1)]..., 1];
+ #r_new = [1, r...]
+ #V = [rand(Uniform(-1, 1), (r_new[i], n[i], r_new[i+1])) for i=1:d];
+ #D = ttmps_eval(V, n)
+ #C, singular_val, errors= tt_svd(D, n, r, d);
+ #p = plot(errors,
+ # ylabel="error",
+ # label=L"$\frac{\|\hat{A}_k - A\|_f}{\|A\|_f}$",
+ # xlabel=L"k",
+ # lw = 1,
+ # xticks=collect(1:d),
+ # markershape=:circle,
+ # legend=:topleft,
+ # margin=5mm,
+ # dpi=300,
+ # title="TT-SVD of uniform tensor")
+ #savefig(p, "./plots/ttsvd-uniform-error.png")
+
+ ## exercise 7
- # Exercise 3
d = 4; n = 51;
+ ndim = [n for i=1:d]
A = [f([t1, t2, t3, t4]) for t1 in t(n), t2 in t(n), t3 in t(n), t4 in t(n)];
B = [g([t1, t2, t3, t4]) for t1 in t(n), t2 in t(n), t3 in t(n), t4 in t(n)];
- # Exercise 3a
- σ_unf_a = σ_unfold(A, d);
- σ_unf_b = σ_unfold(B, d);
- p1 = plot(margin=5mm)
- p2 = plot(margin=5mm)
- for k=1:d
- plot!(p1,
- σ_unf_a[k],
- ylabel=L"$\log(\Sigma_{A_k})$",
- label=L"$\Sigma_{A_%$k}$",
- xlabel=L"n",
- lw = 1,
- yaxis=:log,
- markershape=:circle,
- dpi=300,
- title="Unfolding of A")
- plot!(p2,
- σ_unf_b[k],
- ylabel=L"$\log(\Sigma_{B_k})$",
- label=L"$\Sigma_{_B%$k}$",
- xlabel=L"n",
- lw = 1,
- dpi=300,
- markershape=:circle,
- yaxis=:log,
- title="Unfolding of B")
- end
- savefig(p1, "./plots/singular-dec-A.png")
- savefig(p2, "./plots/singular-dec-B.png")
+ ## Exercise 7a
- ϵ_s = [1/(10^(j*2)) for j=1:5] # computer not goode enough for j = 6
-
- r_jk_a, ϵ_jk_a, σ_jk_a = rank_approx(A);
- r_jk_b, ϵ_jk_b, σ_jk_b = rank_approx(B);
+ r_jk_a, ϵ_jk_a, σ_jk_a = rank_approx(A, ttmps_unfold);
+ r_jk_b, ϵ_jk_b, σ_jk_b = rank_approx(B, ttmps_unfold);
p1 = plot(size=[700, 350], margin=5mm, dpi=300)
p2 = plot(size=[700, 350], margin=5mm, dpi=300)
for k=1:d
@@ -102,122 +249,8 @@ function main()
legend=:topleft,
title=L"Tensor $g(x_1, x_2, x_3, x_4)$")
end
- savefig(p1, "./plots/hosvd-error-A.png")
- savefig(p2, "./plots/hosvd-error-B.png")
-
- # Exercise 3b
- Nj_hosvd_a = []; ϵ_hosvd_a = []; σ_hosvd_a = [];
- Nj_hosvd_b = []; ϵ_hosvd_b = []; σ_hosvd_b = [];
-
- println("Checking validity of error analysis for HOSVD")
- for j=1:length(ϵ_s)
- S0_a, Vs_a, vals_a, errors_a = my_hosvd(A, r_jk_a[j, :])
- N_vals = 0
- for V in Vs_a
- N_vals += length(V)
- end
- append!(Nj_hosvd_a, N_vals + length(S0_a))
- append!(ϵ_hosvd_a, errors_a[end])
- append!(σ_hosvd_a, [vals_a])
- println("ϵ_hosvd_a[j] <= |ϵ_jk_a|_f : $(ϵ_hosvd_a[j] <= norm(ϵ_jk_a)) ϵ_hosvd_a[j] = $(ϵ_hosvd_a[j])")
-
- S0_b, Vs_b, vals_b, errors_b = my_hosvd(B, r_jk_b[j, :])
- N_vals = 0
- for V in Vs_b
- N_vals += length(V)
- end
- append!(Nj_hosvd_b, N_vals + length(S0_b))
- append!(ϵ_hosvd_b, errors_b[end])
- append!(σ_hosvd_b, [vals_b])
- println("ϵ_hosvd_b[j] <= |ϵ_jk_b|_f : $(ϵ_hosvd_b[j] <= norm(ϵ_jk_b)) ϵ_hosvd_b[j] = $(ϵ_hosvd_b[j])")
- end
-
- # exercise 3b
- j = 5
- p1 = plot(dpi=300, size=[800, 350], margin=5mm)
- p2 = plot(dpi=300, size=[800, 350], margin=5mm)
- for k=1:d
- plot!(p1, σ_unf_a[k][1:r_jk_a[j, k]] ./ σ_hosvd_a[j][k],
- ylabel=L"\frac{\sigma^{A_k}_\alpha}{\sigma^{hosvd}_{k\alpha}}",
- label=L"k=%$k",
- xlabel=L"j",
- xticks=collect(1:length(r_jk_a)),
- yticks=[1, 1.001],
- ylim=[1, 1.001],
- lw = 2,
- markershape=:circle,
- legend=:topleft,
- title=L"$f(x_1, x_2, x_3, x_4)$")
- plot!(p2, σ_unf_b[k][1:r_jk_b[j, k]] ./ σ_hosvd_b[j][k],
- ylabel=L"\frac{\sigma^{B_k}_\alpha}{\sigma^{hosvd}_{k\alpha}}",
- label=L"k=%$k",
- xlabel=L"j",
- yticks=[1, 1.001],
- ylim=[1, 1.001],
- lw = 2,
- markershape=:circle,
- legend=:topleft,
- title=L"$g(x_1, x_2, x_3, x_4)$")
- end
- savefig(p1, "./plots/hosvd-sigmaratio-a.png")
- savefig(p2, "./plots/hosvd-sigmaratio-b.png")
-
- # Exercise 3d
-
- p1 = plot(size=[700, 350], margin=5mm, dpi=300)
- p2 = plot(size=[700, 350], margin=5mm, dpi=300)
- plot!(p1,
- 1 ./ ϵ_s,
- Nj_hosvd_a,
- ylabel=L"$N_{j}$",
- label=L"A-$N_{j}$",
- xlabel=L"\varepsilon_{j}^{-1}",
- lw = 2,
- markershape=:circle,
- xaxis=:log10,
- legend=:topleft,
- title=L"Tucker of $f(x_1, x_2, x_3, x_4)$")
- plot!(p2,
- 1 ./ ϵ_s,
- Nj_hosvd_b,
- ylabel=L"$N_{j}$",
- label=L"B-$N_{j}$",
- xlabel=L"\varepsilon_{j}^{-1}",
- lw = 2,
- markershape=:circle,
- xaxis=:log10,
- legend=:topleft,
- title=L"Tucker of $g(x_1, x_2, x_3, x_4)$")
- savefig(p1, "./plots/hosvd-Nj-a.png")
- savefig(p2, "./plots/hosvd-Nj-b.png")
-
- # exercise 6
- d = 4;
- n = [20+k for k=1:d];
- r = [[2*k for k=1:(d-1)]..., 1];
- r_new = [1, r...]
- V = [rand(Uniform(-1, 1), (r_new[i], n[i], r_new[i+1])) for i=1:d];
- D = ttmps_eval(V, n)
- C, singular_val, errors= tt_svd(D, n, r, d);
- p = plot(errors,
- ylabel="error",
- label=L"$\frac{\|\hat{A}_k - A\|_f}{\|A\|_f}$",
- xlabel=L"k",
- lw = 1,
- xticks=collect(1:d),
- markershape=:circle,
- legend=:topleft,
- margin=5mm,
- dpi=300,
- title="TT-SVD of uniform tensor")
- savefig(p, "./plots/ttsvd-uniform-error.png")
-
- # exercise 7
-
- d = 4; n = 51;
- ndim = [n for i=1:d]
- A = [f([t1, t2, t3, t4]) for t1 in t(n), t2 in t(n), t3 in t(n), t4 in t(n)];
- B = [g([t1, t2, t3, t4]) for t1 in t(n), t2 in t(n), t3 in t(n), t4 in t(n)];
+ savefig(p1, "./plots/tucker-error-A.png")
+ savefig(p2, "./plots/tucker-error-B.png")
# exercise 7b
Nj_ttsvd_a = []; ϵ_ttsvd_a = []; σ_ttsvd_a = [];
@@ -275,6 +308,7 @@ function main()
savefig(p1, "./plots/ttsvd-sigmaratio-a.png")
savefig(p2, "./plots/ttsvd-sigmaratio-b.png")
+
# Exercise 7d
p1 = plot(size=[700, 350], margin=5mm, dpi=300)
p2 = plot(size=[700, 350], margin=5mm, dpi=300)
diff --git a/sesh4/src/plots/hosvd-Nj-b.png b/sesh4/src/plots/hosvd-Nj-b.png
Binary files differ.
diff --git a/sesh4/src/plots/hosvd-error-B.png b/sesh4/src/plots/hosvd-error-B.png
Binary files differ.
diff --git a/sesh4/src/plots/hosvd-sigmaratio-a.png b/sesh4/src/plots/hosvd-sigmaratio-a.png
Binary files differ.
diff --git a/sesh4/src/plots/hosvd-sigmaratio-b.png b/sesh4/src/plots/hosvd-sigmaratio-b.png
Binary files differ.
diff --git a/sesh4/src/plots/hosvd-uniform-error.png b/sesh4/src/plots/hosvd-uniform-error.png
Binary files differ.
diff --git a/sesh4/src/plots/ttsvd-uniform-error.png b/sesh4/src/plots/ttsvd-uniform-error.png
Binary files differ.
diff --git a/sesh4/src/unf-aprx.jl b/sesh4/src/unf-aprx.jl
@@ -18,12 +18,29 @@ function σ_unfold(C, d)
return Σ_s
end
+function ttmps_unfold(A, k)
+ n = size(A)
+ d = length(n)
+ C = copy(A)
+ A_k = reshape(C, (prod([n[i] for i in 1:k]), prod([n[i] for i in k+1:d])))
+ return A_k
+end
+
+function tucker_unfold(A, k)
+ n = size(A)
+ d = length(n)
+ C = copy(A)
+ C_perm = permutedims(C, ([i for i=1:d if i!=k]..., k))
+ A_k = reshape(C_perm, (prod([n[i] for i=1:d if i!=k]), n[k]))
+ return A_k
+end
+
ϵ_s = [1/(10^(j*2)) for j=1:5] # computer not goode enough for j = 6
-function rank_approx(C, ϵ_s=ϵ_s)
+function rank_approx(C, method,ϵ_s=ϵ_s)
d = length(size(C))
ϵ_jk = []; r_jk = []; σ_jk = []
for k=1:d
- C_k = tenmat(C, k)
+ C_k = method(C, k)
for (j, ϵ_j) in enumerate(ϵ_s)
for r=1:rank(C_k)
U_hat, Σ_hat, V_hat = tsvd(C_k, r)
diff --git a/sesh5/src/arithmetic.jl b/sesh5/src/arithmetic.jl
@@ -0,0 +1,93 @@
+#/usr/bin/julia
+
+using LinearAlgebra
+using Plots
+using Plots.Measures
+using Distributions
+using LaTeXStrings
+using TSVD: tsvd
+using TensorToolbox: tenmat, ttm, contract # todo:implement
+
+@doc raw"""
+ Addition of two tensors in the TT-MPS format of $A$ and $B$ of the same size
+ $n_1 × ⋯ × n_d$ with different TT-MPS decompositions U with ranks $p_1, … , p_{d-1}$
+ and V with ranks $q_1, … ,q_{d-1}$ respectively.
+
+ Output is the TT-MPS representation of $A+B$ with ranks $(q_1 + p_1), ⋯ , (q_{d-1} + p_{d-1})$.
+"""
+function tt_add(U, V, α, β)
+ W = []
+ d = length(U)
+ for k ∈ 1:d
+ U_k = U[k]; V_k = V[k]
+ p_k, n_k, p_k1 = size(U_k)
+ q_k, n_k, q_k1 = size(V_k)
+ if k == 1
+ W_1 = zeros(1, n_k, q_k1 + p_k1)
+ for i=1:n_k
+ W_1[:, i, :] += hcat(α * U_k[:, i, :], β * V_k[:, i, :])
+ end
+ append!(W, [W_1])
+ elseif k == d
+ W_d = zeros(q_k + p_k, n_k, 1)
+ for i=1:n_k
+ W_d[:, i, :] += vcat(U_k[:, i, :], V_k[:, i, :])
+ end
+ append!(W, [W_d])
+ else
+ W_k = zeros(p_k + q_k, n_k, q_k1 + p_k1)
+ for i=1:n_k
+ W_k[:, i, :] += hcat(vcat(U_k[:, i, :], zeros(q_k, p_k1)),
+ vcat(zeros(p_k, q_k1), V_k[:, i, :]))
+ end
+ append!(W, [W_k])
+ end
+ end
+ return W
+end
+
+@doc raw"""
+ Hadamard product of two tensors in the TT-MPS format of $A$ and $B$ of the same size
+ $n_1 × ⋯ × n_d$ with different TT-MPS decompositions U with ranks $p_1, … , p_{d-1}$
+ and V with ranks $q_1, … ,q_{d-1}$ respectively.
+
+ Output is the representation of $A ⊙ B$ with ranks $(q_1 ⋅ p_1), ⋯ , (q_{d-1} ⋅ p_{d-1})$.
+"""
+function tt_mult(U, V)
+ W = []
+ d = length(U)
+ for k ∈ 1:d
+ p_k, n_k, p_k1 = size(U[k])
+ q_k, n_k, q_k1 = size(V[k])
+ W_k = zeros(p_k*q_k, n_k, p_k1 * q_k1)
+ for i_k ∈ 1:n_k
+ W_k[:, i_k, :] += kron(U[k][:, i_k, :], V[k][:, i_k, :])
+ end
+ append!(W, [W_k])
+ end
+ return W
+end
+
+@doc raw"""
+ Matrix Vector product of two tensors in the TT-MPS format of $A$, size $(m_1⋯m_d×n_1⋯n_d)$ and $u$ size $(n_1⋯n_d)$
+ with different TT-MPS decompositions U with ranks $p_1, … , p_{d-1}$
+ and V with ranks $q_1, … ,q_{d-1}$ respectively.
+
+ Output is the TT-MPS representation of $Au$ with ranks $(q_1 ⋅ p_1), ⋯ , (q_{d-1} ⋅ p_{d-1})$.
+"""
+function tt_matvec(A, U)
+ W = []
+ d = length(U)
+ for k ∈ 1:d
+ α_k, m_k, n_k, α_k1 = size(A[k])
+ β_k, m_k, β_k1 = size(U[k])
+ W_k = zeros(α_k * β_k, n_k, α_k1 * β_k1)
+ for i_k ∈ 1:m_k
+ for j_k ∈ 1:n_k
+ W_k[:, i_k, :] += A[k][:, i_k, j_k, :] * U[k][:, j_k, :]
+ end
+ end
+ append!(W, [W_k])
+ end
+ return W
+end
diff --git a/sesh5/src/helpers.jl b/sesh5/src/helpers.jl
@@ -0,0 +1,26 @@
+#/usr/bin/julia
+
+using LinearAlgebra
+using Plots
+using Plots.Measures
+using Distributions
+using LaTeXStrings
+using TSVD: tsvd
+using TensorToolbox: tenmat, ttm, contract # todo:implement
+
+@doc raw"""
+ Helper functions to calculated the tensor values of the Chebyshev Polynomials \textit{che_pol}
+ on a Grid \textit{t(n)}
+"""
+
+t(n) = [2*(i-1)/(n-1)-1 for i ∈ 1:n]
+
+function che_pol(x, q)
+ if abs(x) ≤ 1
+ return cos(q*acos(x))
+ elseif x ≥ 1
+ return cosh(q*acosh(x))
+ elseif x ≤ -1
+ return (-1)^q * cosh(q*acosh(-x))
+ end
+end
diff --git a/sesh5/src/main.jl b/sesh5/src/main.jl
@@ -0,0 +1,91 @@
+#/usr/bin/julia
+
+using LinearAlgebra
+using Plots
+using Plots.Measures
+using Distributions
+using LaTeXStrings
+using TSVD: tsvd
+using TensorToolbox: tenmat, ttm, contract # todo:implement
+
+include("arithmetic.jl")
+include("helpers.jl")
+include("ttmps.jl")
+
+function main()
+ # Exercise 3a)
+ n = 10
+ d = 4
+ p_Σ = []
+ p_R = []
+ for (p, q) ∈ ([3, 4], [5, 7])
+ X = [che_pol(sum([t_1/1, t_2/2, t_3/3, t_4/4]), p) for t_1 ∈ t(n), t_2 ∈ t(n), t_3 ∈ t(n), t_4 ∈ t(n)];
+ Y = [che_pol(sum([t_1/1, t_2/2, t_3/3, t_4/4]), q) for t_1 ∈ t(n), t_2 ∈ t(n), t_3 ∈ t(n), t_4 ∈ t(n)];
+ S = X .+ Y;
+ Z = X .* Y;
+ global X_Σ = []; global Y_Σ = []; global S_Σ = []; global Z_Σ = []
+ global r_x = []; global r_y = []; global r_s = []; global r_z = []
+ for k ∈ 1:d-1
+ X_k = ttmps_unfold(X, k)
+ Y_k = ttmps_unfold(Y, k)
+ S_k = ttmps_unfold(S, k)
+ Z_k = ttmps_unfold(Z, k)
+ _, σ_x, _ = svd(X_k)
+ _, σ_y, _ = svd(Y_k)
+ _, σ_s, _ = svd(S_k)
+ _, σ_z, _ = svd(Z_k)
+ append!(X_Σ, [σ_x])
+ append!(Y_Σ, [σ_y])
+ append!(S_Σ, [σ_s])
+ append!(Z_Σ, [σ_z])
+ append!(r_x, [rank(X_k)])
+ append!(r_y, [rank(Y_k)])
+ append!(r_s, [rank(S_k)])
+ append!(r_z, [rank(Z_k)])
+ end
+ p_σ = plot(layout = 4, margin=5mm, size=(800, 800), title="p=$p, q=$q", dpi=500)
+ plot!(p_σ[1], X_Σ, markershape=:circle, ylabel=L"\Sigma_{X_k}", label=[L"X_1" L"X_2" L"X_3"])
+ plot!(p_σ[2], Y_Σ, markershape=:circle, ylabel=L"\Sigma_{Y_k}", label=[L"Y_1" L"Y_2" L"Y_3"])
+ plot!(p_σ[3], S_Σ, markershape=:circle, ylabel=L"\Sigma_{S_k}", label=[L"S_1" L"S_2" L"S_3"])
+ plot!(p_σ[4], Z_Σ, markershape=:circle, ylabel=L"\Sigma_{Z_k}", label=[L"Z_1" L"Z_2" L"Z_3"])
+ append!(p_Σ, [p_σ])
+ p_r = plot(margin=5mm, size=(700, 350), title="p=$p, q=$q", ylabel=L"r_k", dpi=300)
+ plot!(p_r, r_x, markershape=:circle, label=L"r_X")
+ plot!(p_r, r_y, markershape=:circle, label=L"r_Y")
+ plot!(p_r, r_s, markershape=:circle, label=L"r_S", alpha=0.5)
+ plot!(p_r, r_z, markershape=:circle, label=L"r_Z")
+ append!(p_R, [p_r])
+
+ end
+ savefig(p_Σ[1], "./plots/sigma_34.png")
+ savefig(p_Σ[2], "./plots/sigma_57.png")
+ savefig(p_R[1], "./plots/rank_34.png")
+ savefig(p_R[2], "./plots/rank_57.png")
+
+
+ # Exercise 3b)
+ n = 10; p = 5; q = 7
+ dims = [10 for i ∈ 1:4]
+
+ X = [che_pol(sum([t_1/1, t_2/2, t_3/3, t_4/4]), p) for t_1 ∈ t(n), t_2 ∈ t(n), t_3 ∈ t(n), t_4 ∈ t(n)];
+ Y = [che_pol(sum([t_1/1, t_2/2, t_3/3, t_4/4]), q) for t_1 ∈ t(n), t_2 ∈ t(n), t_3 ∈ t(n), t_4 ∈ t(n)];
+ S = X .+ Y
+ Z = X .* Y
+
+ C_x, r_x, σ_x, ϵ_x = tt_svd_ϵ(X, 1e-12);
+ C_y, r_y, σ_y, ϵ_y = tt_svd_ϵ(Y, 1e-12);
+ C_s = tt_add(C_x, C_y, 1, 1);
+ C_z = tt_mult(C_x, C_y);
+
+ Q_s, r_s = t_mpstt_ϵ(C_s, 1e-12);
+ Q_z, r_z = t_mpstt_ϵ(C_z, 1e-7);
+
+ println("Exercise 3b")
+ println("|ttmps_eval(C_s) - S||_F/||S||_F = $(norm(S - ttmps_eval(C_s, dims))/norm(S))")
+ println("|ttmps_eval(C_z) - Z||_F/||Z||_F = $(norm(Z - ttmps_eval(C_z, dims))/norm(Z))")
+ println()
+ println("|ttmps_eval(Q_s) - S||_F/||S||_F = $(norm(S - ttmps_eval(Q_s, dims))/norm(S))")
+ println("|ttmps_eval(Q_z) - Z||_F/||Z||_F = $(norm(Z - ttmps_eval(Q_z, dims))/norm(Z))")
+end
+
+main()
diff --git a/sesh5/src/plots/rank_34.png b/sesh5/src/plots/rank_34.png
Binary files differ.
diff --git a/sesh5/src/plots/rank_57.png b/sesh5/src/plots/rank_57.png
Binary files differ.
diff --git a/sesh5/src/plots/sigma_34.png b/sesh5/src/plots/sigma_34.png
Binary files differ.
diff --git a/sesh5/src/plots/sigma_57.png b/sesh5/src/plots/sigma_57.png
Binary files differ.
diff --git a/sesh5/src/ttmps.jl b/sesh5/src/ttmps.jl
@@ -0,0 +1,152 @@
+#/usr/bin/julia
+
+using LinearAlgebra
+using Plots
+using Plots.Measures
+using Distributions
+using LaTeXStrings
+using TSVD: tsvd
+using TensorToolbox: tenmat, ttm, contract # todo:implement
+
+@doc raw"""
+ Output: the k-th TT-MPS unfolding matrix of $A$.
+ A_{i_1, … , i_n} -> A_{i_1 ⋯ i_{k}, i_{k+1} ⋯ i_d}
+"""
+function ttmps_unfold(A, k)
+ n = size(A)
+ d = length(n)
+ C = copy(A)
+ A_k = reshape(C, (prod([n[i] for i in 1:k]), prod([n[i] for i in k+1:d])))
+ return A_k
+end
+
+@doc raw"""
+ TT-MPS decomposition evaluation to a Tensor
+"""
+function ttmps_eval(U, n)
+ A = U[1]
+ for U_k ∈ U[2:end]
+ A = contract(A, U_k)
+ end
+ return reshape(A, n...)
+end
+
+@doc raw"""
+ Truncated MPS-TT:
+ Given a decomposition U of ranks $p_1, … ,p_{d-1}$ produce a decomposition of target
+ ranks not exceeding $r = [r_1, … ,r_{d-1}]$ with the TT-MPS orthogonalization algorithm
+"""
+function t_mpstt(U, r)
+ r_n = [1, r..., 1]
+ d = length(U)
+ Q = []; U_k = U[1]
+ for k ∈ 2:d
+ α_k_1, i_k, α_k = size(U_k)
+ U_k_bar = reshape(U_k, (α_k_1*i_k, α_k))
+
+ P_k, Σ_k, W_k = tsvd(U_k_bar, r_n[k])
+ Q_k = reshape(P_k, (α_k_1, i_k, r_n[k]))
+ Z_k = Diagonal(Σ_k) * transpose(W_k)
+ U_k = contract(Z_k, U[k])
+ append!(Q, [Q_k])
+ if k == d
+ append!(Q, [U_k])
+ end
+ end
+ return Q
+end
+
+@doc raw"""
+ Truncated MPS-TT:
+ Same as above but with a tolerace = ϵ, with which the ranks are
+ approximated. The difference is there are no target ranks required
+"""
+function t_mpstt_ϵ(U, ϵ)
+ r = [1]
+ d = length(U)
+ dims = [size(U[k], 2) for k ∈ 1:d]
+ δ = ϵ/sqrt(d-1) * norm(ttmps_eval(U, dims))
+ Q = []; U_k = U[1]
+ for k ∈ 2:d
+ α_k_1, i_k, α_k = size(U_k)
+ U_k_bar = reshape(U_k, (α_k_1*i_k, α_k))
+ for r_k ∈ 1:rank(U_k_bar)
+ global P_k, Σ_k, W_k = tsvd(U_k_bar, r_k)
+ U_k_bar_hat = P_k * Diagonal(Σ_k) * transpose(W_k)
+ if norm(U_k_bar - U_k_bar_hat)/norm(U_k_bar) ≤ δ
+ append!(r, r_k)
+ break
+ end
+ end
+ Q_k = reshape(P_k, (α_k_1, i_k, r[k]))
+ Z_k = Diagonal(Σ_k) * transpose(W_k)
+ U_k = contract(Z_k, U[k])
+ append!(Q, [Q_k])
+ if k == d
+ append!(Q, [U_k])
+ end
+ end
+ return Q, r
+end
+
+
+@doc raw"""
+ TT-SVD algorithm
+"""
+function tt_svd(A, r)
+ n = size(A)
+ d = length(n)
+ r_new = [1, r...]
+ S_0_hat = copy(A)
+ S0s = []
+ C = []; σ = []; ϵ = []
+ for k=2:d
+ B_k = reshape(S_0_hat, (r_new[k-1] * n[k-1], prod([n[i] for i=k:d])))
+ U_hat, Σ_hat, V_hat = tsvd(convert(Matrix{Float64}, B_k), r_new[k])
+ C_k = reshape(U_hat, (r_new[k-1], n[k-1], r_new[k]))
+ W_k_hat = Diagonal(Σ_hat) * transpose(V_hat)
+ S_0_hat = reshape(W_k_hat, (r_new[k], [n[i] for i=k:d]...))
+
+ append!(C, [C_k])
+ append!(σ, [Σ_hat])
+ A_hat = ttmps_eval([C..., S_0_hat], n)
+ append!(ϵ, norm(A_hat - A)/norm(A))
+ end
+ append!(C, [reshape(S_0_hat, (size(S_0_hat)..., 1))])
+ return C, σ, ϵ
+end
+
+@doc raw"""
+ Tolerace bound TT-SVD algorithm, no target ranks required.
+ The ranks are calculated based on the tolerace specified.
+"""
+function tt_svd_ϵ(A, tol)
+ n = size(A)
+ d = length(n)
+ r = [1]
+ S_0_hat = copy(A)
+ δ = tol/sqrt(d-1) * norm(A)
+ C = []; σ = []; ϵ = []
+ for k=2:d
+ B_k = reshape(S_0_hat, (r[k-1] * n[k-1], prod([n[i] for i=k:d])))
+ for r_k = 1:rank(B_k)
+ global U_hat, Σ_hat, V_hat = tsvd(convert(Matrix{Float64}, B_k), r_k)
+ B_k_hat = U_hat * Diagonal(Σ_hat) * transpose(V_hat)
+ if norm(B_k - B_k_hat)/norm(B_k) ≤ δ
+ append!(r, r_k)
+ break
+ end
+ end
+ C_k = reshape(U_hat, (r[k-1], n[k-1], r[k]))
+ W_k_hat = Diagonal(Σ_hat) * transpose(V_hat)
+ S_0_hat = reshape(W_k_hat, (r[k], [n[i] for i=k:d]...))
+
+ append!(C, [C_k])
+ append!(σ, [Σ_hat])
+ A_hat = ttmps_eval([C..., S_0_hat], n)
+ append!(ϵ, norm(A_hat - A)/norm(A))
+ end
+ append!(C, [reshape(S_0_hat, (size(S_0_hat)..., 1))])
+ return C, r, σ, ϵ
+end
+
