Contents

Generate Data Set

clear;
% dimension, np, nq
m = 500, np = 1500, nq = 1500
rng(10); adj = rand(m)<.001; adj = triu(adj, 1) + triu(adj, 1)'; thetaP = double(adj)*.45 + eye(m)*2;
rng(10); adj = rand(m)<.001-3e-5; adj = triu(adj, 1) + triu(adj, 1)'; thetaQ = double(adj)*.45 + eye(m)*2;
rng(1);
xp = mvnrnd(zeros(m,1),inv(thetaP),np)'; xq = mvnrnd(zeros(m,1),inv(thetaQ),nq)';
kp = kernel_linear(xp); kq = kernel_linear(xq);
m =

   500


np =

        1500


nq =

        1500

Naive subgradient descent

theta = sparse(zeros(size(kq,1),1)); lambda = 0.55*log(m)/sqrt(np);
tic
step = 1; slength = inf; iter = 0; fold = inf;
while(slength > 1e-5)
    [f, gt] = LLKLIEP(theta,kp,kq);
    g = zeros(size(gt));

    id = abs(theta)>0;
    g(id) = gt(id) + lambda*sign(theta(id));
    id = theta==0 & gt > lambda;
    g(id) = gt(id) - lambda;
    id = theta==0 & gt < -lambda;
    g(id) = gt(id) + lambda;
    theta = theta - step*g./(iter+1);
    slength = step*norm(g)./(iter+1);
    fdiff = abs(f - fold);

    %display some stuffs
    if iter > 5000
        disp('max iteration reached.')
        break;
    else
        iter = iter+1;
        fdiff = abs(f - fold);
        fold = f;
        if ~mod(iter,100)
            disp(sprintf('%d, %.5f, %.5f, %.5f, nz: %d',...
                iter, slength,fdiff,full(fold),sum(theta(1:end-m)~=0)))
        end
    end
end
toc
100, 0.00025, 0.00005, 7.23998, nz: 5
200, 0.00011, 0.00002, 7.23684, nz: 5
300, 0.00006, 0.00001, 7.23528, nz: 5
400, 0.00004, 0.00001, 7.23428, nz: 5
500, 0.00003, 0.00001, 7.23356, nz: 5
600, 0.00003, 0.00000, 7.23300, nz: 5
700, 0.00002, 0.00000, 7.23256, nz: 5
800, 0.00002, 0.00000, 7.23218, nz: 5
900, 0.00002, 0.00000, 7.23187, nz: 5
1000, 0.00001, 0.00000, 7.23159, nz: 5
1100, 0.00001, 0.00000, 7.23135, nz: 5
1200, 0.00001, 0.00000, 7.23113, nz: 5
1300, 0.00001, 0.00000, 7.23094, nz: 5
时间已过 805.260965 秒。

visualize

tt = ones(m); tt = triu(tt,1); tt=tt(:); idx = tt~=0;
Delta = zeros(1,m*m);
Delta(idx) = theta(1:end-m); Delta = reshape(Delta,m,m);
Delta = Delta + Delta';

hh = figure; spy(thetaP - thetaQ, 64, 'r'); hold on;
spy(abs(Delta), 18); title(sprintf('lambda = %.5f',lambda))
h = legend('ground truth', 'detected'); h.Location='southeast';

Function LLKLIEP

function [l,g,h] = LLKLIEP(theta,kP,kQ)
l = -mean(theta'*kP,2) + log(sum(exp(theta'*kQ),2));
N_q = sum(exp(theta'*kQ),2);
g_q = exp(theta'*kQ)./ N_q;
g = -mean(kP,2) + kQ*g_q';
% hessian
if nargout>2
    HH = diag(g_q) - g_q'*g_q;
    h = kQ*HH*kQ';
end
end

Function kernel_linear

function [k] = kernel_linear(x)
[d n] = size(x); k = [];
tt = ones(d); tt = triu(tt,1); tt=tt(:); idx = tt~=0;
for i = 1:n
    t = x(:,i) * x(:,i)';
    t = triu(t,1); t = t(:);
    k = [k, t(idx)];
end
k = [k;x.^2];
end