ivan@152
|
1 function DL=SMALL_two_step_DL(Problem, DL)
|
ivan@152
|
2
|
ivan@152
|
3 % determine which solver is used for sparse representation %
|
ivan@152
|
4
|
ivan@152
|
5 solver = DL.param.solver;
|
ivan@152
|
6
|
bmailhe@209
|
7 % determine which type of udate to use ('KSVD', 'MOD', 'ols', 'opt' or 'LGD') %
|
ivan@152
|
8
|
ivan@152
|
9 typeUpdate = DL.name;
|
ivan@152
|
10
|
ivan@152
|
11 sig = Problem.b;
|
ivan@152
|
12
|
ivan@152
|
13 % determine dictionary size %
|
ivan@152
|
14
|
ivan@152
|
15 if (isfield(DL.param,'initdict'))
|
ivan@152
|
16 if (any(size(DL.param.initdict)==1) && all(iswhole(DL.param.initdict(:))))
|
ivan@152
|
17 dictsize = length(DL.param.initdict);
|
ivan@152
|
18 else
|
ivan@152
|
19 dictsize = size(DL.param.initdict,2);
|
ivan@152
|
20 end
|
ivan@152
|
21 end
|
ivan@152
|
22 if (isfield(DL.param,'dictsize')) % this superceedes the size determined by initdict
|
ivan@152
|
23 dictsize = DL.param.dictsize;
|
ivan@152
|
24 end
|
ivan@152
|
25
|
ivan@152
|
26 if (size(sig,2) < dictsize)
|
ivan@152
|
27 error('Number of training signals is smaller than number of atoms to train');
|
ivan@152
|
28 end
|
ivan@152
|
29
|
ivan@152
|
30
|
ivan@152
|
31 % initialize the dictionary %
|
ivan@152
|
32
|
luis@198
|
33 if (isfield(DL.param,'initdict'))
|
ivan@152
|
34 if (any(size(DL.param.initdict)==1) && all(iswhole(DL.param.initdict(:))))
|
ivan@152
|
35 dico = sig(:,DL.param.initdict(1:dictsize));
|
ivan@152
|
36 else
|
ivan@152
|
37 if (size(DL.param.initdict,1)~=size(sig,1) || size(DL.param.initdict,2)<dictsize)
|
ivan@152
|
38 error('Invalid initial dictionary');
|
ivan@152
|
39 end
|
ivan@152
|
40 dico = DL.param.initdict(:,1:dictsize);
|
ivan@152
|
41 end
|
ivan@152
|
42 else
|
ivan@152
|
43 data_ids = find(colnorms_squared(sig) > 1e-6); % ensure no zero data elements are chosen
|
ivan@152
|
44 perm = randperm(length(data_ids));
|
ivan@152
|
45 dico = sig(:,data_ids(perm(1:dictsize)));
|
ivan@152
|
46 end
|
ivan@152
|
47
|
ivan@152
|
48 % flow: 'sequential' or 'parallel'. If sequential, the residual is updated
|
ivan@152
|
49 % after each atom update. If parallel, the residual is only updated once
|
ivan@152
|
50 % the whole dictionary has been computed. Sequential works better, there
|
ivan@152
|
51 % may be no need to implement parallel. Not used with MOD.
|
ivan@152
|
52
|
ivan@152
|
53 if isfield(DL.param,'flow')
|
ivan@152
|
54 flow = DL.param.flow;
|
ivan@152
|
55 else
|
ivan@152
|
56 flow = 'sequential';
|
ivan@152
|
57 end
|
ivan@152
|
58
|
ivan@152
|
59 % learningRate. If the type is 'ols', it is the descent step of
|
ivan@152
|
60 % the gradient (typical choice: 0.1). If the type is 'mailhe', the
|
bmailhe@209
|
61 % descent step is the optimal step*rho (typical choice: 1, although 2 works
|
bmailhe@209
|
62 % better). Not used for MOD and KSVD.
|
ivan@152
|
63
|
ivan@152
|
64 if isfield(DL.param,'learningRate')
|
ivan@152
|
65 learningRate = DL.param.learningRate;
|
ivan@152
|
66 else
|
bmailhe@209
|
67 switch typeUpdate
|
bmailhe@209
|
68 case 'ols'
|
bmailhe@209
|
69 learningRate = 0.1;
|
bmailhe@209
|
70 otherwise
|
bmailhe@209
|
71 learningRate = 1;
|
bmailhe@209
|
72 end
|
ivan@152
|
73 end
|
ivan@152
|
74
|
ivan@152
|
75 % number of iterations (default is 40) %
|
ivan@152
|
76
|
ivan@152
|
77 if isfield(DL.param,'iternum')
|
ivan@152
|
78 iternum = DL.param.iternum;
|
ivan@152
|
79 else
|
ivan@152
|
80 iternum = 40;
|
ivan@152
|
81 end
|
ivan@152
|
82 % determine if we should do decorrelation in every iteration %
|
ivan@152
|
83
|
ivan@152
|
84 if isfield(DL.param,'coherence')
|
ivan@152
|
85 decorrelate = 1;
|
ivan@152
|
86 mu = DL.param.coherence;
|
ivan@152
|
87 else
|
ivan@152
|
88 decorrelate = 0;
|
ivan@152
|
89 end
|
ivan@152
|
90
|
ivan@152
|
91 % show dictonary every specified number of iterations
|
ivan@152
|
92
|
ivan@153
|
93 if isfield(DL.param,'show_dict')
|
ivan@152
|
94 show_dictionary=1;
|
ivan@152
|
95 show_iter=DL.param.show_dict;
|
ivan@152
|
96 else
|
ivan@152
|
97 show_dictionary=0;
|
ivan@152
|
98 show_iter=0;
|
ivan@152
|
99 end
|
ivan@152
|
100
|
ivan@152
|
101 % This is a small patch that needs to be resolved in dictionary learning we
|
ivan@152
|
102 % want sparse representation of training set, and in Problem.b1 in this
|
ivan@152
|
103 % version of software we store the signal that needs to be represented
|
ivan@152
|
104 % (for example the whole image)
|
ivan@152
|
105
|
ivan@152
|
106 tmpTraining = Problem.b1;
|
ivan@152
|
107 Problem.b1 = sig;
|
ivan@153
|
108 if isfield(Problem,'reconstruct')
|
ivan@153
|
109 Problem = rmfield(Problem, 'reconstruct');
|
ivan@153
|
110 end
|
ivan@152
|
111 solver.profile = 0;
|
ivan@152
|
112
|
ivan@152
|
113 % main loop %
|
ivan@152
|
114
|
ivan@152
|
115 for i = 1:iternum
|
ivan@153
|
116 Problem.A = dico;
|
ivan@152
|
117 solver = SMALL_solve(Problem, solver);
|
luis@198
|
118 [dico, solver.solution] = dico_update(dico, sig, solver.solution, ...
|
luis@198
|
119 typeUpdate, flow, learningRate);
|
luis@198
|
120 if (decorrelate)
|
bmailhe@209
|
121 dico = dico_decorr_symetric(dico, mu, solver.solution);
|
luis@198
|
122 end
|
ivan@153
|
123
|
ivan@152
|
124 if ((show_dictionary)&&(mod(i,show_iter)==0))
|
ivan@152
|
125 dictimg = SMALL_showdict(dico,[8 8],...
|
ivan@152
|
126 round(sqrt(size(dico,2))),round(sqrt(size(dico,2))),'lines','highcontrast');
|
ivan@152
|
127 figure(2); imagesc(dictimg);colormap(gray);axis off; axis image;
|
ivan@152
|
128 pause(0.02);
|
ivan@152
|
129 end
|
ivan@152
|
130 end
|
ivan@152
|
131
|
ivan@152
|
132 Problem.b1 = tmpTraining;
|
ivan@152
|
133 DL.D = dico;
|
ivan@152
|
134
|
ivan@153
|
135 end
|
ivan@153
|
136
|
ivan@153
|
137 function Y = colnorms_squared(X)
|
ivan@153
|
138
|
ivan@153
|
139 % compute in blocks to conserve memory
|
ivan@153
|
140 Y = zeros(1,size(X,2));
|
ivan@153
|
141 blocksize = 2000;
|
ivan@153
|
142 for i = 1:blocksize:size(X,2)
|
ivan@153
|
143 blockids = i : min(i+blocksize-1,size(X,2));
|
ivan@153
|
144 Y(blockids) = sum(X(:,blockids).^2);
|
ivan@153
|
145 end
|
ivan@153
|
146
|
luis@198
|
147 end |