minimax_sq_c_exp_anis_blocks.m

function [lambdas,omegas,Kt,Sm] = minimax_sq_c_exp_anis_blocks(sigma, block_sizes, mu, x, y, iters, tol_lambda, tol_opt, tol_Q, tol_fun, sigma0, ...
    lambda0, dc, dc_iters, dc_draws)

% Best convex combination over a continuous set of anisotropic Gaussian kernels for square loss
% through iterative optimization wrt c and lambda
% sigma : #blocks x 2 matrix of [sigma_min, sigma_max]
% sigma0 : n x #blocks matrix
% block_sizes : sizes for each block of consecutive vector components. Should sum to d.
% tol_lambda, tol_opt, tol_Q, tol_fun : fractional tolerances for Newton search, dc optimization, 
%                                     quadratic form and functional respectively

if (nargin == 12)
    dc = 0;
end

n = length(lambda0);
[m,d] = size(x);
B = length(block_sizes);
if (abs(sum(lambda0)-1) > eps)
    disp('minimax_sq_c_exp_anis_blocks() : \t Lambdas should sum to 1');
    return;
end
if (sum(block_sizes)~=d)
    disp('minimax_sq_c_exp_anis_blocks() : \t Block sizes should sum to d');
    return;
end
lambdas = lambda0;
omega0 = 1./sigma0.^2;
omegas = omega0';
omega_min = 1./sigma(:,1).^2;
omega_max = 1./sigma(:,2).^2;

% Compute matrices (x_ik-x_jk)^2
block_starts = [1,1+cumsum(block_sizes)];
for k=1:B
    Dx(k,:,:) = dist_matrix(x(:,block_starts(k):block_starts(k+1)-1),x(:,block_starts(k):block_starts(k+1)-1));
end
Kt = zeros(m,m);
for i=1:n
    Kt = Kt + lambda0(i)*exp(-squeeze(sum(Dx.*repmat(omega0(i,:)',[1,m,m]),1)));
end
Sm = inf;

% Nested function definitions for DC programming
    function gval = g_quad(w)
        gval = (exp(-w'*D1))*C1;
    end
    function hval = h_quad(w)
        hval = (exp(-w'*D2))*C2;
    end
    function aval = a_feas(w)
        aval = max(max([omega_min-w,w-omega_max]));
    end
    function ggval = gradg(w)
        ggval = -(repmat(C1',B,1).*D1)*(exp(-D1'*w));
    end
    function ghval = gradh(w)
        ghval = -(repmat(C2',B,1).*D2)*(exp(-D2'*w));
    end
    function gaval = suba(w)
        [m_min,imin] = max(omega_min-w);
        [m_max,imax] = max(w-omega_max);
        gaval = zeros(B,1);
        if (m_min > m_max)
            gaval(imin) = -1;
        else
            gaval(imax) = 1;
        end
    end
    function gminval = gmin(gfn,S,tol)
        VP = S{1};
        [dS,vp] = size(VP);
        cP = S{2};
        x0 = VP(:,floor(rand(1)*vp)+1);
        options = optimset('TolFun',tol*abs(Q),'TolX',tol*norm(x0),'TolCon',tol*norm(x0),'Display','off',...
	  'LargeScale','off');
        [wopt,gminval] = fmincon(gfn,x0,cP(1:dS,:)',-cP(dS+1,:)',[],[],[],[],[],options);
    end




    function [hval,ghval] = ht(w)
        hval = -(exp(-w(1:B)'*D2))*C2+w(B+1);
        ghval = [(repmat(C2',B,1).*D2)*(exp(-D2'*w(1:B)));1];
    end
    function [c,ceq,gc,gceq] = con(w)
        c = [];
        gc = [];
        ceq = (exp(-w(1:B)'*D1))*C1-w(B+1);
        gceq = [-(repmat(C1',B,1).*D1)*(exp(-D1'*w(1:B)));-1];
    end




for t=1:iters
    ch = -2*mu*inv(Kt+mu*eye(m))*y;
    Q = ch'*Kt*ch;

    % Maximization of (c,G_x(omega)c)

    if (dc)
        Ch = -ch*ch';
        ind1 = find(Ch > 0);
        ind2 = find(Ch < 0);
        l1 = length(ind1);
        l2 = length(ind2);
        C1 = Ch(ind1);
        C2 = -Ch(ind2);
        D1 = Dx(:,ind1);
        D2 = Dx(:,ind2);
        % The initial simplex is generated by the lower bounds hyperplanes
        % and the hyperplane 1'*w - 1'*omega_max = 0
        smax = sum(omega_max);
        S0{1} = [omega_min,repmat(omega_min,1,B)+(smax-sum(omega_min))*eye(B)];
        S0{2} = [-eye(B),ones(B,1);omega_min',-smax];
        S0{3} = [ones(1,B),0;ones(B,B)-eye(B),ones(B,1)];
        S0{4} = ones(B+1,B+1)-eye(B+1);
        gmax = inline('max(gfn(S{1}))','gfn','S','tol');
        if (t==1)
            w0 = (omega_min+omega_max)/2;
        else
            w0 = omegas(:,1);
        end
        [opt,minval] = cutting_plane_fmincon_vertex(@g_quad,@h_quad,@a_feas,w0,@gradg,@gradh,@suba,S0,@gmin,gmax,...
            tol_opt,tol_lambda,dc_iters,dc_draws);

if 0
%         t0 = feval(@g_quad,w0);
%         options2 = optimset('TolFun',tol_opt*abs(Q),'TolX',tol_opt*norm([w0;t0]),'TolCon',tol_opt*norm([w0;t0]),...
%             'Display','off','LargeScale','off','GradConstr','on','GradObj','on');
%         [opt,minval,exitflag] = fmincon(@ht,[w0;t0],[],[],[],[],[],[],@con,options2);
% 	opt = opt(1:B);
% 	fprintf('exitflag = %d \n',exitflag);
end



        fmax = -minval;
        fprintf('opt = ');
        for r=1:B
            fprintf('%.8g ',opt(r));
        end
        fprintf('\t fmax = %.8g \t',fmax);
        fprintf('Q = %.8g \n',Q);

        if 0
        % DEBUG
        intervals = 30;
        [pts1,pts2] = meshgrid(exp(log(omega_min(1)):(log(omega_max(1))-log(omega_min(1)))/intervals:log(omega_max(1))),...
            exp(log(omega_min(2)):(log(omega_max(2))-log(omega_min(2)))/intervals:log(omega_max(2))));
        for p1=1:intervals+1
            for p2=1:intervals+1
                res(p1,p2)=g_quad([pts1(p1,p2);pts2(p1,p2)])-h_quad([pts1(p1,p2);pts2(p1,p2)]);
            end
        end
        [mm,i1] = min(res);
        [mm,i2] = min(mm);
        figure; hold on;
        contour(pts1,pts2,res);
        plot(pts1(i1(i2),i2),pts2(i1(i2),i2),'x');
        fprintf('min value = %.8g \n',res(i1(i2),i2));
        [cntr,hcntr] = contour(pts1,pts2,res,[Q Q]);
        clabel(cntr,hcntr,[Q]);
        set(gca,'XScale','log');
        set(gca,'YScale','log');
        % END DEBUG
        end
        
    else
        % ********* UNTESTED ***********
        %         intervals_min = omega_min;
        %         intervals_max = omega_max;
        %         init_range = omega_max-omega_min;
        %         fmax = -inf;
        %         while (intervals_min)
        %             a = intervals_min(:,1);
        %             b = intervals_max(:,1);
        %             [im,in] = size(intervals_min);
        %             if (max((b-a)./init_range) < tol_omega)
        %                 intervals_min = intervals_min(:,[2:in]);
        %                 intervals_max = intervals_max(:,[2:in]);
        %                 continue;
        %             end
        %             mid = (a+b)/2;
        %             % Gradient ascent
        %             omegan = mid;
        %             covered = [mid,mid];
        %             temp = exp(-dist_matrix2(x,x,diag(mid)));
        %             fn = ch'*temp*ch;
        %             % Gradient computation
        %             for k=1:d
        %                 dfn(k,1) = -ch'*(squeeze(Dx(k,:,:)).*temp)*ch;
        %             end
        %             for k=1:d
        %                 for l=k:d
        %                     J(k,l) = ch'*(squeeze(Dx(k,:,:)).*squeeze(Dx(l,:,:)).*temp)*ch;
        %                     J(l,k) = J(k,l);
        %                 end
        %             end
        %             if (cond(J) > 1/eps)
        %                 intervals_min = [a,mid,intervals_min(:,[2:in])];
        %                 intervals_max = [mid,b,intervals_max(:,[2:in])];
        %                 if (fn > fmax)
        %                     opt = mid;
        %                     fmax = fn;
        %                 end
        %                 continue;
        %             end
        %             eta = abs((inv(J)*dfn)./dfn);
        %             while (max(abs(dfn.*eta)./init_range) < tol_omega)
        %                 eta = 2*eta;
        %             end
        %             step = dfn.*eta;
        %             while (max(abs(step./init_range)) > tol_omega)
        %                 if (find(omegan+step < a | omegan+step > b))
        %                     eta = eta/2;
        %                     if (max(eta) < eps | max(abs(step/2)./init_range) < tol_omega)
        %                         left = find(omegan+step < a);
        %                         covered(left,1) = a(left);
        %                         right = find(omegan+step > b);
        %                         covered(right,2) = b(right);
        %                         break;
        %                     end
        %                 else
        %                     temp = ch'*exp(-dist_matrix2(x,x,diag(omegan+step)))*ch;
        %                     if (temp < fn)
        %                         eta = eta/2;
        %                         if (max(eta) < eps | max(abs(step/2)./init_range) < tol_omega)
        %                             break;
        %                         end
        %                     else
        %                         omegan = omegan + step;
        %                         fn = temp;
        %                         covered(:,1) = min([covered(:,1),omegan],[],2);
        %                         covered(:,2) = max([covered(:,2),omegan],[],2);
        %                         % DEBUG
        %                         plot(omegan(1),omegan(2),'ro');
        %                         % END DEBUG
        %                     end
        %                 end
        %                 for k=1:d
        %                     dfn(k) = -ch'*(squeeze(Dx(k,:,:)).*exp(-dist_matrix2(x,x,diag(omegan))))*ch;
        %                 end
        %                 if (max(abs(dfn)) < eps)
        %                     covered(:,1) = max([min([covered(:,1),omegan-init_range*tol_omega],[],2),a],[],2);
        %                     covered(:,2) = min([max([covered(:,2),omegan+init_range*tol_omega],[],2),b],[],2);
        %                     break;
        %                 end
        %                 step = dfn.*eta;
        %             end
        %             % Compare with value in memory
        %             if (fn > fmax)
        %                 opt = omegan;
        %                 fmax = fn;
        %                 if (fmax > g*(1+tol_e))
        %                     break;
        %                 end
        %             end
        %             % Exclude covered interval
        %             next = [a,b];
        %             for k=1:d
        %                 temp1 = next;
        %                 temp2 = next;
        %                 temp1(k,2) = covered(k,1);
        %                 temp2(k,1) = covered(k,2);
        %                 intervals_min = [temp1(:,1),temp2(:,1),intervals_min(:,[2:in])];
        %                 intervals_max = [temp1(:,2),temp2(:,2),intervals_max(:,[2:in])];
        %                 next(k,:) = covered(k,:);
        %             end
        %             % DEBUG
        %             line(covered(:,1),[covered(1,1),covered(2,2)],'k');
        %             line([covered(1,1),covered(2,2)],covered(:,2),'k');
        %             line(covered(:,2),[covered(1,2),covered(2,1)],'k');
        %             line([covered(1,2),covered(2,1)],covered(:,1),'k');
        %             %plot(covered,[fn,fn],'kx');
        %             % END DEBUG
        %         end
    end

    % DEBUG
    %plot(opt(1),opt(2),'ro');
    % END DEBUG

    if (fmax <= Q*(1+tol_Q))       % Optimum reached (no kernel can improve)
        fprintf('minimax_sq_c_exp_anis_blocks() : \t converged in %d iterations\n\n',t);
        return;
    end

    %DEBUG
    %disp(sprintf('opt = %.3g \n',opt));

    Kj = exp(-squeeze(sum(Dx.*repmat(opt,[1,m,m]),1)));

    % DEBUG
    %     figure;
    %     hold on;
    %     for l=0:0.01:1
    %         plot(l,-mu*y'*inv(l*Kj+(1-l)*Kt+mu*eye(m))*y);
    %     end
    % END DEBUG

    % Search for lambda minimizing phi
    l_prev = inf;
    l = 0.5;
    while (abs(l-l_prev) > tol_lambda)
        tempi = inv(l*Kj+(1-l)*Kt+mu*eye(m));
        cKl = -2*mu*tempi*y;
        temp = (Kt-Kj)*cKl;
        dphi = cKl'*temp/(4*mu);
        ddphi = temp'*tempi*temp/(2*mu);
        l_prev = l;
        l = l-dphi/ddphi;
        if (l>=1)
            l=1;
        elseif (l<=0)
            l=0;
        end
    end

    %plot(l,-mu*y'*inv(l*Kj+(1-l)*Kt+mu*eye(m))*y,'rx');
    % update kernel
    Kt = l*Kj+(1-l)*Kt;
    omegas = [opt,omegas];
    lambdas = [l,(1-l)*lambdas];

    %disp(sprintf('l = %f \n',l));
    Sm_prev = Sm;
    Sm = mu*y'*inv(Kt+mu*eye(m))*y;
    if (abs(Sm-Sm_prev)/Sm < tol_fun)
        fprintf('minimax_sq_c_exp_anis_blocks() : \t converged in %d iterations\n',t);
        fprintf('minimax_sq_c_exp_anis_blocks() : \t fmax-Q = %.3g \t Sm = %.3g  dSm = %.3g \n\n',fmax-Q,Sm,Sm_prev-Sm);
        return;
    end
    fprintf('minimax_sq_c_exp_anis_blocks() : \t iteration %d \t fmax-Q = %.3g \t Sm = %.3g  dSm = %.3g \n',t,fmax-Q,Sm,Sm_prev-Sm);
end

fprintf('minimax_sq_c_exp_anis_blocks() : \t could not converge after %d iterations\n\n',t);
end