scm_infsup_lb.m

Go to the documentation of this file.

function alpha_LB = scm_infsup_lb(mu,model,operator,SCM_data)
%alpha_LB = scm_infsup_lb(mu,model,model_data)
% computation of the lowerbound of the coercive constant 
% equal to the coercive constant computet via SCM
%
% mu         : parameter which lower bound value is desired
% model      : standart model 
% model_data : standart model data
% SCM_data   : precomputed offline SCM data (e.g. via scm_coercive_greedy)

% David Kreplin 01.12 2014


if isfield(SCM_data,'C') == 0
    error('SCM Greedy computation is required first')
end

ind =@ (i,j,q) j + (i-1)*q - i*(i-1)/2;
[~,Ind_PMa]   = scm_get_neighbours(SCM_data.params.Ma, mu, SCM_data.C.Set);
[PMp,Ind_PMp] = scm_get_neighbours(SCM_data.params.Mp, mu, SCM_data.params.ParameterSet);
Ma = length(Ind_PMa);
Mp = length(Ind_PMp);
b_lp= zeros(Ma+Mp,1);
A_lp= zeros(Ma+Mp,(SCM_data.params.Q)*(SCM_data.params.Q+1)/2);

b_lp(1:Ma) = -1*SCM_data.C.infsup(Ind_PMa(1:Ma));
A_lp(1:Ma,:) = SCM_data.C.ParameterVector(Ind_PMa(:,1:Ma),:);

for n = 1:Mp
    model = set_mu(model,PMp(:,n));
    model.decomp_mode=2;
    [A,~] =  operator(model);
    
    for i=1:SCM_data.params.Q
        for j=i:SCM_data.params.Q
            A_lp(n+Ma,ind(i,j,model.SCM_data.Q)) = -(2-(i==j))*(A(i)*A(j));
        end
    end
end

model = set_mu(model,mu);
model.decomp_mode=2;
[A,~] =  operator(model);
f_lp = zeros(1,(SCM_data.params.Q)*(SCM_data.params.Q+1)/2);
for i=1:SCM_data.params.Q
    for j=i:SCM_data.params.Q
        f_lp(ind(i,j,SCM_data.params.Q)) = (2-(i==j))*(A(i)*A(j));
    end
end

Options = optimset('Display','off');
%keyboard;
[~,alpha_LB]=linprog(f_lp,A_lp,b_lp,[],[],SCM_data.lower_bound_infsup,...
    SCM_data.upper_bound_infsup,SCM_data.upper_bound_infsup,Options);
%keyboard;
end