2 %model=scm_greedy(model, model_data)
5 % computes necessary values
for the calculation of the infsup constant
6 % requires model and model data structure
12 % SCM_data.Ma : number of parameters with an already
13 % computed constant near to mu
14 % SCM_data.Mp : number of parameters without computed
16 % SCM_data.NumberOfParameters : size of the test set used in the SCM
18 % SCM_data.eps : error boundary of the
Greedy procedure
19 % SCM_data.max_iterations : maximal number of
Greedy iterations
20 % SCM_data.parallel : = 0 non parallel version
21 % = 1 parallel (pool stays active after
23 % = 2 parallel (closing pool afterwards)
24 % David Kreplin 05.01.2015
27 warning('No SCM properties are set -> using standart values')
31 params.NumberOfParameters = 100;
33 params.max_iterations=inf;
37 if ~isfield(params, 'Ma')
41 if ~isfield(params, 'Mp')
45 if ~isfield(params, 'NumberOfParameters')
46 warning('SCM: No Number of Parameters is set -> setting NumberOfParameters = 100');
47 params.NumberOfParameters = 100;
50 if ~isfield(params, 'eps')
54 if ~isfield(params, 'max_iterations')
55 params.max_iterations = inf;
58 if ~isfield(params, 'parallel')
62 if ~isfield(params, 'para_rand_seed')
63 params.para_rand_seed = 1234;
68 if ~isfield(params, 'ParameterSet')
69 rand('state',params.para_rand_seed);
70 params.ParameterSet = rand_uniform(params.NumberOfParameters,model.mu_ranges);
72 NumberOfParameters = length(params.ParameterSet(1,:));
74 SCM_data.params = params;
76 ind =@ (i,j,q) j + (i-1)*q - i*(i-1)/2;
81 %% Boundary BOX for YLB
83 [A,~] = operator(model);
84 A=remove_dirichlet(A,1,SCM_data.params.dirichlet_gids );
88 SCM_data.lower_bound_infsup = zeros(Q*(Q+1)/2,1);
89 SCM_data.upper_bound_infsup = zeros(Q*(Q+1)/2,1);
91 disp('-- SCM: computing boundary box');
93 %innerproduct_matrix = model.get_inner_product_matrix(model_data);
94 innerproduct_matrix(SCM_data.params.dirichlet_gids ,:)=[];
95 innerproduct_matrix(:,SCM_data.params.dirichlet_gids )=[];
96 innerproduct_matrix = 0.5*(innerproduct_matrix + innerproduct_matrix');
100 fun = @(x) A{1,i}
'*(innerproduct_matrix\(A{1,j}*x));
102 [lambda1, success] = safeEig(fun, size(A{i}, 1), innerproduct_matrix, 'lr
');
104 [lambda2, success] = safeEig(fun, size(A{i}, 1), innerproduct_matrix, 'sr
');
109 spectrum = eig(full(A{1,i}'*(innerproduct_matrix
'\A{1,j})), full(innerproduct_matrix));
110 if any(abs(imag(spectrum)) > 1e-9)
113 spectrum = real(spectrum);
114 lambda1 = max(spectrum);
115 lambda2 = min(spectrum);
118 SCM_data.upper_bound_infsup(ind(i,j,Q)) = lambda1;
119 SCM_data.lower_bound_infsup(ind(i,j,Q)) = lambda2;
123 % SCM_data.lower_bound_infsup=[0;0;-7.3747e-16];
124 % SCM_data.upper_bound_infsup=[1;0;1.267];
129 SCM_data.C=AddC(params.ParameterSet(:,1),model,operator,innerproduct_matrix,SCM_data);
130 SCM_data.params.ParameterSet(:,1)=[];
131 NumberOfParameters = NumberOfParameters - 1;
136 %% SCM Greedy - parallel version
137 if SCM_data.params.parallel >= 1
138 % s = matlabpool('size
'); %OLDER MATLAB VERSION
140 % matlabpool('open
');
142 if isempty(gcp('nocreate
'))
145 %reduce parallel data amount
146 ParameterSet = SCM_data.params.ParameterSet;
147 while SCM_data.params.eps <= eta
149 eta_temp=zeros(NumberOfParameters,1);
150 parfor i = 1:NumberOfParameters
151 alpha_UB = scm_infsup_ub(ParameterSet(:,i),model,operator,SCM_data);
152 alpha_LB = scm_infsup_lb(ParameterSet(:,i),model,operator,SCM_data);
153 eta_temp(i) = (sqrt(alpha_UB)-sqrt(alpha_LB))/sqrt(alpha_UB);
156 [eta,index]=max(eta_temp);
158 SCM_data.C=AddC(ParameterSet(:,index),model,operator,innerproduct_matrix,SCM_data);
159 ParameterSet(:,index) = [];
160 NumberOfParameters = NumberOfParameters-1;
162 disp(['-- SCM: step
' num2str(n) ' - error:
' num2str(eta)]);
163 if n >= SCM_data.params.max_iterations
167 SCM_data.params.ParameterSet=ParameterSet;
168 if params.parallel == 2
169 % matlabpool('close
'); %OLDER MATLAB VERSION
170 delete(gcp('nocreate
'))
174 %% non-parallel version
175 while SCM_data.params.eps <= eta
177 eta_temp=zeros(NumberOfParameters,1);
178 for i = 1:NumberOfParameters
179 alpha_UB = scm_infsup_ub(SCM_data.params.ParameterSet(:,i),model,operator,SCM_data);
180 alpha_LB = scm_infsup_lb(SCM_data.params.ParameterSet(:,i),model,operator,SCM_data);
181 eta_temp(i) = abs(sqrt(alpha_UB)-sqrt(alpha_LB))/sqrt(alpha_UB);
183 [eta,index]=max(eta_temp);
185 SCM_data.C=AddC(SCM_data.params.ParameterSet(:,index),model,operator,innerproduct_matrix,SCM_data);
187 SCM_data.params.ParameterSet(:,index) = [];
188 NumberOfParameters = NumberOfParameters-1;
190 disp(['-- SCM: step
' num2str(n) ' - error:
' num2str(eta)]);
191 if n >= SCM_data.params.max_iterations
203 function [lambda, success] = safeEig(fun, sizeA, B, sigma)
209 while ~success && opts.tol < 1e-2
212 lambda = eigs(fun, sizeA, B, 1, sigma, opts);
215 warning('safeEig: eigs exception, lowering the tolerance
');
216 opts.tol = 100*opts.tol;
220 if abs(imag(lambda)) > 1e-9
221 warning('safeEig: ignoring complex eigenvalues
');
223 lambda = real(lambda);
228 function C=AddC(mu,model,operator,innerproduct_matrix,SCM_data)
229 ind =@ (i,j,q) j + (i-1)*q - i*(i-1)/2;
230 % innerproduct_matrix = model.get_inner_product_matrix(model_data);
231 % innerproduct_matrix(SCM_data.params.dirichlet_gids ,:)=[];
232 % innerproduct_matrix(:,SCM_data.params.dirichlet_gids )=[];
234 if isempty(SCM_data.C)
236 model = set_mu(model,mu);
238 [A,~] = operator(model);
239 A=remove_dirichlet(A,0,SCM_data.params.dirichlet_gids);
240 %[V,D] = eig(full(A')*inv(full(innerproduct_matrix))
'*full(A),full(innerproduct_matrix));
241 %[v,d] = fast_eigs(A'*(innerproduct_matrix
'\A),innerproduct_matrix,'sm
' );
243 fun = @(X) A
'\(innerproduct_matrix*(A\X));
244 [v,d] = eigs(fun, size(A, 1), innerproduct_matrix, 1, 'sm
');
246 C.infsup =d; %infsup;
250 [A_Matrix,~] = operator(model);
251 A_Matrix=remove_dirichlet(A_Matrix,1,SCM_data.params.dirichlet_gids );
252 y_UB=zeros(SCM_data.params.Q*(SCM_data.params.Q+1)/2,1);
253 norm_value = v'*innerproduct_matrix*v;
254 for i=1:SCM_data.params.Q
255 for j=i:SCM_data.params.Q
256 y_UB(ind(i,j,SCM_data.params.Q)) = v
'*A_Matrix{1,i}'*(innerproduct_matrix
'\(A_Matrix{1,j}*v))/norm_value;
262 [A,~] = operator(model);
263 ParameterVector=zeros(0,(SCM_data.params.Q)*(SCM_data.params.Q+1)/2);
264 for i=1:SCM_data.params.Q
265 for j=i:SCM_data.params.Q
266 ParameterVector(ind(i,j,SCM_data.params.Q)) = -(2-(i==j))*(A(i)*A(j));
269 C.ParameterVector = ParameterVector;
271 %Parameter for lower bound computation
275 C.Set = [SCM_data.C.Set mu];
276 model = set_mu(model,mu);
278 [A,~] = operator(model);
279 A=remove_dirichlet(A,0,SCM_data.params.dirichlet_gids);
280 % [V,D] = eig(full(A')*inv(full(innerproduct_matrix))
'*full(A),full(innerproduct_matrix));
281 %[v,d] = fast_eigs(A'*(innerproduct_matrix
'\A),innerproduct_matrix,'sm
' );
283 fun = @(X) A
'\(innerproduct_matrix*(A\X));
284 [v,d] = eigs(fun, size(A, 1), innerproduct_matrix, 1, 'sm
');
286 %C.infsup = [C.infsup;infsup];
287 C.infsup = [SCM_data.C.infsup;d];
291 [A_Matrix,~] = operator(model);
292 A_Matrix=remove_dirichlet(A_Matrix,1,SCM_data.params.dirichlet_gids);
293 y_UB=zeros(SCM_data.params.Q*(SCM_data.params.Q+1)/2,1);
294 norm_value = v'*innerproduct_matrix*v;
295 for i=1:SCM_data.params.Q
296 for j=i:SCM_data.params.Q
297 y_UB(ind(i,j,SCM_data.params.Q)) = v
'*A_Matrix{1,i}'*(innerproduct_matrix
'\(A_Matrix{1,j}*v))/norm_value;
302 [A,~] = operator(model);
303 ParameterVector=zeros(0,(SCM_data.params.Q)*(SCM_data.params.Q+1)/2);
304 for i=1:SCM_data.params.Q
305 for j=i:SCM_data.params.Q
306 ParameterVector(ind(i,j,SCM_data.params.Q)) = -(2-(i==j))*(A(i)*A(j));
309 C.ParameterVector = [SCM_data.C.ParameterVector;ParameterVector];
311 C.UBVector = [SCM_data.C.UBVector y_UB];
317 %% Upperbound computation
319 function alpha_UB = scm_infsup_ub(mu,model,operator,SCM_data)
320 ind =@ (i,j,q) j + (i-1)*q - i*(i-1)/2;
322 model = set_mu(model,mu);
324 [A,~] = operator(model);
325 factor= zeros(1,(SCM_data.params.Q)*(SCM_data.params.Q+1)/2);
327 for i=1:SCM_data.params.Q
328 for j=i:SCM_data.params.Q
329 factor(ind(i,j,SCM_data.params.Q)) = (2-(i==j))*(A(i)*A(j));
333 alpha_UB=min(factor*SCM_data.C.UBVector);
341 function A=remove_dirichlet(A,mode,I)
Customizable implementation of an abstract greedy algorithm.
function SCM_data = scm_infsup_greedy(model, params, operator, innerproduct_matrix)
model=scm_greedy(model, model_data)