-
Notifications
You must be signed in to change notification settings - Fork 0
/
epsilon_runs_rand_range.m
176 lines (138 loc) · 6.08 KB
/
epsilon_runs_rand_range.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
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
function [results,info] = epsilon_runs_rand_range(X,A,krng,srng,crng,trng,ntrials)
%EPSILON_RUNS Run experiments on SAMPLED_SOLVE
%
% RESULTS = EPSILON_RUNS(X,A,KRNG,SRNG,CRNG,TRNG,NT) runs a series of
% experiments on tensor X with the factor matrices in cell array A.
% - KRNG = array of modes to use
% - SRNG = array of number of samples (s) to use
% - CRNG = cell array of functions to compute c as a function of s,
% where c is the max expected copies of any sample
% - TRNG = cell array of functions to compute t as a function of s,
% where t is the threshold for deterministic samples
% - NT = number of trials for each combo
%
% RESULTS = EPSILON_RUNS(X,R,KRNG,...) differs in that R is an
% integer rank for a random initial guess.
%%
fprintf('\n*** Epsilon Runs ***\n');
%% Set up for the tensor
fprintf('\n---Extracting data from tensor---\n');
% Size
d = ndims(X);
sz = size(X);
fprintf('Tensor is of size %s with nnz = %d\n', tt_intvec2str(sz), nnz(X));
%%
fprintf('\n---Calculating fiber indices---\n');
tic
Xfidxs = zeros(nnz(X),d);
for kk = 1:d
Xfidxs(:,kk) = fiber_indices(X,kk);
end
toc
%% Set up solution
fprintf('\n---Setup Factor Matrices---\n')
if isscalar(A)
r = A;
A = cell(d, 1);
fprintf('\tUsing randomized range finder\n');
rnginit;
for kk = 1:d
A{kk} = randomized_range_finder(X, kk, Xfidxs(:, kk), r, 100000);
end
else
fprintf('Using provided factor matrices\n');
r = size(A{1},2);
end
fprintf('Rank = %d\n', r);
%% Leverage scores
fprintf('\n---Calculating leverage scores---\n');
tic
Alev = cell(d,1);
for kk = 1:d
Alev{kk} = tt_leverage_scores(A{kk});
end
toc
%% Setup data
nk = length(krng); % Number of k-values
ns = length(srng); % Number of s-values
nc = length(crng); % Number of c-functions
nt = length(trng); % Number of t-functions
full_solve_time = zeros(nk,1);
results = cell(nk, ns, nc, nt, ntrials);
%% Mode outer loop
for kidx = 1:nk
k = krng(kidx);
%% Full solve for mode k
fprintf('\n---Full Solve for k=%d---\n',k);
tic
Vf = full_solve(X,A,k); %<-- This would be way faster with MTTKRP pre-computed
full_solve_time(kidx) = toc;
%% Sample loop
for sidx = 1:ns
s = srng(sidx);
for cidx = 1:nc
cfunc = crng{cidx};
c = cfunc(s);
for tidx = 1:nt
tfunc = trng{tidx};
t = tfunc(s);
for trialidx = 1:ntrials
fprintf('\n---Sampled Solve for k=%d, s=%d, c=%g, t=%g, trial=%d---\n', k, s, c, t, trialidx);
rnginit;
% Compute sampled solve
tic
[Vs,info] = sampled_solve(X,A,Alev,k,s,c,t,Xfidxs(:,k),true,true);
itr.sampled_solve_time_plus_post = toc;
% Extract sample-sample residual
itr.ss_resid = info.ss_residual;
% Extract timings
itr.timings = info.timings;
itr.sampled_solve_time = info.totaltime;
itr.posttime = info.posttime;
% Extract some info about the damping, weights, and repeats
itr.maxwgt = max(info.wgts);
itr.minwgt = min(info.wgts);
% Extra damping factors
itr.alpha = info.alpha;
% Extract most frequent sampled indices
[midx_uniq,~,mm] = unique(info.midx,'row');
midx_cnts = accumarray(mm,1);
midx_aug = [midx_cnts midx_uniq];
midx_foo = sortrows(midx_aug,'descend');
itr.midx_top10 = midx_foo(1:10,:);
% Extract information about deterministic indices
itr.sdet = info.sdet;
itr.pdet = info.pdet;
midx_det_aug = [info.pvec info.midx_det];
midx_det_foo = sortrows(midx_det_aug,'descend');
itr.midx_det_top10 = midx_det_foo(1:min(10, itr.sdet), :);
% Calculate errors
fprintf('Calculating Gamma...\n');
tic
[itr.gamma, itr.ff_resid, itr.fs_resid] = full_relerr(X,k,A,Vf,Vs); %<-- This would be way faster with MTTKRP pre-computed
toc
itr.epsilon = abs(itr.ff_resid - itr.ss_resid)/max(1, itr.ff_resid);
% Print info
fprintf('\tSampled solve time: %fs\n', itr.sampled_solve_time);
fprintf('\tFull-Full residual : %f\n', itr.ff_resid);
fprintf('\tFull-Sparse residual : %f\n', itr.fs_resid);
fprintf('\tSparse-Sparse residual: %f\n', itr.ss_resid);
fprintf('\tGamma : %e\n', itr.gamma);
fprintf('\tEpsilon: %e\n', itr.epsilon);
fprintf('\tDampling Factors: [ ');
fprintf('%f ', itr.alpha);
fprintf(']\n');
fprintf('\tMin/Max weight: %g / %g\n', itr.minwgt, itr.maxwgt);
fprintf('\tTop-10 Most Frequent Samples:\n');
fprintf('\t\t cnt=%d midx=[%d %d %d]\n', itr.midx_top10')
fprintf('\t %i Deterministic Indices with Total Probability %f \n', itr.sdet, itr.pdet);
if itr.sdet > 0
fprintf('\t\t prob=%.3e midx=[%d %d %d]\n', itr.midx_det_top10')
end
% Save everything
results{kidx,sidx,cidx,tidx,trialidx} = itr;
end %trialidx
end % tidx
end % cidx
end % sidx
end % kidx