compute_graph.m

function [con_graph, alpha, A] = compute_graph(W, knn_idx, self_b, cross_b,b)
% Learn a consistent graph from multiple graphs.
% Inputs:
%   W - weight matrix of a graph
%   knn_idx - common kNN index for all views
% Optional Inputs:
%   tol, tol2 - the tolerance that determines convergence of algorithm
% Outputs:
%   con_graph - weight matrix of the learned unified graph
%   E - a cell matrix containing the inconsistent part of all views
%   A - a cell matrix containing the consistent part of all views
% See "Youwei Liang, Dong Huang, and Chang-Dong Wang. Consistency Meets 
% Inconsistency: A Unified Graph Learning Framework for Multi-view 
% Clustering. 2019 IEEE International Conference on Data Mining(ICDM)."
% for a detailed description of the algorithm.
% Author: Youwei Liang
% 2019/08/31

v = length(W);
if nnz(W{1})/numel(W{1}) < 0.4  % if W contains a large proportion of zeros, use sparse mode
    for i=1:v
        W{i} = sparse(W{i});
    end
    sparse_mode = true;
else
    for i=1:v
        W{i} = full(W{i});
    end
    sparse_mode = false;
end
v = length(W);
% b = cross_b*ones(v) - diag(cross_b*ones(1,v)) + diag(self_b*ones(1,v));
b_coef = b + eye(v);
n = size(W{1}, 1);
baW = cell(v,1);
special_baW = cell(v,1);
true_baW = cell(v,1);
A = cell(v,1);
B = cell(v,1);
E = cell(v,1);
up_knn_idx = triu(knn_idx);

zz = 2.^(0:v-1);
ww = 1:2^v-2; % alpha can't be all zeros, so -2
logww = log2(ww);
yy = ww(abs(floor(logww)-logww)>eps);
alpha_zeros_ones = de2bi([0,zz,yy]);
n_eye_coef = -eye(v);

% initialize A{i}, alpha, con_graph
alpha = ones(v,1) / v;
con_graph = W{1};
if sparse_mode
    D = sparse(n, n);
else
    D = zeros(n,n);
end
for i=1:v
    D = max(D, W{i});
    A{i} = full(W{i});
%     temp = A{i}(knn_idx);
%     A{i}(knn_idx) = temp.*(1.5 - rand(m,1));
%     A{i} = (A{i}+A{i}')/2;
%     A{i}(knn_idx) = temp.*(2 - 1.9*rand(m,1));
end
for i=1:v
    if sparse_mode
        A{i} = sparse(A{i});
    end
    A{i} = min(A{i}, D);
end
% con_graph = con_graph/v + rand(m,1);


    % fix A{i}, update con_graph and alpha
    for i=1:v
        E{i} = W{i} - A{i};
    end
    coef = zeros(v);
    coef2 = zeros(v);
    for i=1:v
        for j=i:v
            coef(i,j) = sum(sum(A{i}.*A{j}));
            coef(j,i) = coef(i,j);
            coef2(i,j) = sum(sum(E{i}.*E{j}));
            coef2(j,i) = coef2(i,j);
        end
    end
    coef2 = coef2 .* b;
    
    % compute coefficient for the linear equation
    H = 2*(diag(diag(coef)) - coef/v + coef2);
    one = ones(v, 1);
    
    for i=1:1
        mpl = alpha_zeros_ones(i,:);
        coef3 = H .* ~mpl + n_eye_coef .* mpl;
        X = [coef3, one; 1-mpl, 0];
        temp_b = [zeros(v,1); 1];
        if det(X) == 0 % abs(det(X)) <= eps
%             fprintf('*************')
            solution = pinv(X)*temp_b;
        else
            solution = X \ temp_b;
        end
        alpha = EProjSimplex_new(solution(1:v));
    end
%     fprintf('best_obj:%.3f\n', best_obj)
    con_graph = alpha(1)*A{1};
    for j=2:v
        con_graph = con_graph + alpha(j)*A{j};
    end
    con_graph = con_graph/v;
    
    % fix con_graph and alpha, update A{i}
    alp_coef = alpha * alpha';
    coef = alp_coef .* b_coef;
    if sparse_mode
        commom_baW = sparse(n, n);
    else
        commom_baW = zeros(n,n);
    end
    for i=1:v
        baW{i} = cross_b*alpha(i)*W{i};
        special_baW{i} = self_b*alpha(i)*W{i};
        commom_baW = commom_baW + baW{i};
    end

    for i=1:v
        true_baW{i} = commom_baW-baW{i}+special_baW{i};
        temp = full(alpha(i)*(con_graph + true_baW{i}));
        B{i} = temp(up_knn_idx);
    end
    right_b = cat(2, B{:})';
    if det(coef) == 0
        solution = (pinv(coef) * right_b)';
%         fprintf('------------')
    else
        solution = (coef \ right_b)';
    end
    solution(solution<0) = 0;
    for i=1:v
        temp = solution(:,i);
        A{i} = zeros(n, n);
        A{i}(up_knn_idx) = temp;
        A{i} = max(A{i}, A{i}');
        A{i} = min(W{i}, A{i});
        if sparse_mode
            A{i} = sparse(A{i});
        end
    end

end