forked from kbatseli/TTr1SVD
-
Notifications
You must be signed in to change notification settings - Fork 0
/
getAtilde.m
98 lines (88 loc) · 3.13 KB
/
getAtilde.m
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
function Atilde=getAtilde(U,sigmas,V,sigmaI,n)
% Atilde=getAtilde(U,sigmas,V,sigmaI,n)
% -------------------------------------
% Reconstructs a tensor as a linear combination of rank terms determined
% from a TTr1 decomposition.
%
% Atilde = tensor, the desired tensor from taking a linear combination
% of rank-1 terms,
%
% U = cell, cell of U vectors obtained from the TTr1
% decomposition,
%
% sigmas = vector, contains the singular values of each of the rank-1
% terms in the order as specified in sigmaI,
%
% V = cell, cell of V vectors obtained from the TTr1
% decomposition,
%
% sigmaI = vector, contains indices of rank-1 terms that need to be
% used in the reconstruction. These are the indices of the
% desired leaves of the last level in the TTr1-tree, counting from
% left to right,
%
% n = vector, dimensions of the original tensor A.
%
% Reference
% ---------
%
% A Constructive Algorithm for Decomposing a Tensor into a Finite Sum of Orthonormal Rank-1 Terms
% http://arxiv.org/abs/1407.1593
%
% 2014, Kim Batselier, Haotian Liu, Ngai Wong
% get indices of end node for each level based on size of original tensor n
d=length(n); % order of the tensor
r=ones(1,d); % nodes are paired in groups of r(i) on level i-1
nodesperlevel=ones(1,d); % level i-1 contains nodesperlevel(i) nodes
endI=ones(1,d); % endI(i) is the index of the last SVD of level i-1
for i=2:d
r(i)=min(n(i-1),prod(n(i:end)));
nodesperlevel(i)=prod(r(1:i));
endI(i)=sum(nodesperlevel(1:i));
end
%last endI points to last SVD of previous level
endI=endI(1:end-2);
% initialize output
Atilde=zeros(prod(n),1);
for j=1:length(sigmaI) % for each rank-1 term
tempEnd=endI;
tempR=r;
k=sigmaI(j); % the offset with respect to endI for level d-i, initalize to sigmaI
%% first step is different because we need to multiply with V as well!
% determine new offsets
if mod(k,tempR(end))==0
k=k/tempR(end);
l=tempR(end);
else
l=mod(k,tempR(end));
k=ceil(k/tempR(end));
end
outerprod=kron(V{tempEnd(end)+k}(:,l),U{tempEnd(end)+k}(:,l));
%% now do remaining steps
tempEnd=tempEnd(1:end-1);
tempR=tempR(1:end-1);
for i=1:length(tempEnd)
% determine new offsets
if mod(k,tempR(end))==0
k=k/tempR(end);
l=tempR(end);
else
l=mod(k,tempR(end));
k=ceil(k/tempR(end));
end
outerprod=kron(outerprod,U{tempEnd(end)+k}(:,l));
tempEnd=tempEnd(1:end-1);
tempR=tempR(1:end-1);
end
%% last step is always first node, so k=1
% determine new l-offset
if mod(k,tempR(end))==0
l=tempR(end);
else
l=mod(k,tempR(end));
end
outerprod=kron(outerprod,U{1}(:,l));
Atilde=Atilde+sigmas(j)*outerprod;
end
Atilde=reshape(Atilde,n);
end