链路预测算法MATLAB代码实现与仿真

链路预测，乍听之下有些学术气息，但其核心问题其实非常直观：在一个复杂网络中，哪些尚未连接的节点对之间，未来最有可能产生新的连边？这是复杂网络分析中的经典任务，广泛应用于推荐系统、生物网络分析以及社交网络挖掘等场景。

程序概述

本文提供的MATLAB程序实现了一套完整的链路预测工具箱，覆盖了四大类主流方法：

• 基于局部信息的相似性指标——Common Neighbors、Jaccard、Adamic-Adar 等
• 基于路径的相似性指标——Katz 指数
• 基于随机游走的相似性指标——Rooted PageRank、SimRank
• 矩阵分解方法——潜在特征学习

核心代码实现

下面给出核心类的完整MATLAB实现，涵盖算法调用、性能评估以及可视化展示功能：

classdef LinkPrediction
%LINKPREDICTION 链路预测算法实现
%   包含多种链路预测算法
properties
    A;              % 邻接矩阵
    train_mask;     % 训练掩码矩阵
    test_mask;      % 测试掩码矩阵
    methods;        % 可用的预测方法
end
methods
    function obj = LinkPrediction(adj_matrix, train_ratio)
        %LINKPREDICTION 构造函数
        %   adj_matrix - 网络邻接矩阵
        %   train_ratio - 训练集比例(0-1)
        obj.A = adj_matrix;
        obj.methods = {
            'CommonNeighbors', 'Jaccard', 'AdamicAdar', ...
            'PreferentialAttachment', 'Katz', 'RootedPageRank', ...
            'SimRank', 'MatrixFactorization'
        };
        % 划分训练集和测试集
        if nargin > 1
            obj = obj.split_dataset(train_ratio);
        else
            obj.train_mask = ones(size(obj.A));
            obj.test_mask = zeros(size(obj.A));
        end
    end
    
    function obj = split_dataset(obj, train_ratio)
        %SPLIT_DATASET 划分训练集和测试集
        %   随机隐藏一部分边作为测试集
        [n, ~] = size(obj.A);
        obj.train_mask = ones(n);
        obj.test_mask = zeros(n);
        % 获取所有边的索引
        [rows, cols] = find(triu(obj.A, 1)); % 只取上三角避免重复
        num_edges = length(rows);
        num_train = round(num_edges * train_ratio);
        % 随机选择训练边
        idx = randperm(num_edges);
        train_idx = idx(1:num_train);
        test_idx = idx(num_train+1:end);
        % 创建掩码矩阵
        for i = 1:length(test_idx)
            r = rows(test_idx(i));
            c = cols(test_idx(i));
            obj.train_mask(r, c) = 0;
            obj.train_mask(c, r) = 0;
            obj.test_mask(r, c) = 1;
            obj.test_mask(c, r) = 1;
        end
        % 确保对角线为0
        obj.train_mask(1:n+1:end) = 0;
        obj.test_mask(1:n+1:end) = 0;
    end
    
    function scores = common_neighbors(obj)
        %COMMON_NEIGHBORS 共同邻居算法
        scores = (obj.A * obj.A) .* obj.train_mask;
    end
    
    function scores = jaccard(obj)
        %JACCARD Jaccard相似系数
        [n, ~] = size(obj.A);
        scores = zeros(n);
        for i = 1:n
            for j = i+1:n
                if obj.train_mask(i, j) == 0
                    continue;
                end
                neighbors_i = find(obj.A(i, :));
                neighbors_j = find(obj.A(j, :));
                intersection = length(intersect(neighbors_i, neighbors_j));
                union = length(union(neighbors_i, neighbors_j));
                if union > 0
                    scores(i, j) = intersection / union;
                else
                    scores(i, j) = 0;
                end
                scores(j, i) = scores(i, j);
            end
        end
    end
    
    function scores = adamic_adar(obj)
        %ADAMIC_ADAR Adamic-Adar指数
        [n, ~] = size(obj.A);
        scores = zeros(n);
        % 计算每个节点的度
        degrees = sum(obj.A, 2);
        for i = 1:n
            for j = i+1:n
                if obj.train_mask(i, j) == 0
                    continue;
                end
                common_neighbors = find(obj.A(i, :) & obj.A(j, :));
                score = 0;
                for k = common_neighbors
                    if degrees(k) > 1 % 避免除以0
                        score = score + 1 / log(degrees(k));
                    end
                end
                scores(i, j) = score;
                scores(j, i) = score;
            end
        end
    end
    
    function scores = preferential_attachment(obj)
        %PREFERENTIAL_ATTACHMENT 优先连接算法
        degrees = sum(obj.A, 2);
        scores = (degrees * degrees') .* obj.train_mask;
    end
    
    function scores = katz(obj, beta)
        %KATZ Katz指数
        %   beta - 衰减因子，默认0.01
        if nargin < 2
            beta = 0.01;
        end
        [n, ~] = size(obj.A);
        I = eye(n);
        scores = beta * obj.A; % 长度为1的路径
        % 计算Katz指数：S = βA + β²A² + β³A³ + ...
        % 使用矩阵求逆近似：S = (I - βA)^-1 - I
        scores = inv(I - beta * obj.A) - I;
        scores = scores .* obj.train_mask;
    end
    
    function scores = rooted_pagerank(obj, alpha, max_iter, tol)
        %ROOTED_PAGERANK Rooted PageRank算法
        %   alpha - 随机游走概率，默认0.85
        %   max_iter - 最大迭代次数，默认100
        %   tol - 收敛容差，默认1e-6
        if nargin < 2
            alpha = 0.85;
        end
        if nargin < 3
            max_iter = 100;
        end
        if nargin < 4
            tol = 1e-6;
        end
        [n, ~] = size(obj.A);
        scores = zeros(n);
        % 创建列随机矩阵（转移概率矩阵）
        P = obj.A ./ sum(obj.A, 1);
        P(isnan(P)) = 0; % 处理度为0的节点
        % 对每个节点计算Rooted PageRank
        for i = 1:n
            r = zeros(n, 1);
            r(i) = 1;
            for iter = 1:max_iter
                r_new = alpha * P * r + (1 - alpha) * r;
                if norm(r_new - r, 1) < tol
                    break;
                end
                r = r_new;
            end
            scores(:, i) = r;
        end
        scores = scores .* obj.train_mask;
    end
    
    function scores = simrank(obj, C, max_iter, tol)
        %SIMRANK SimRank算法
        %   C - 衰减因子，默认0.8
        %   max_iter - 最大迭代次数，默认10
        %   tol - 收敛容差，默认1e-4
        if nargin < 2
            C = 0.8;
        end
        if nargin < 3
            max_iter = 10;
        end
        if nargin < 4
            tol = 1e-4;
        end
        [n, ~] = size(obj.A);
        S = eye(n); % 初始化SimRank矩阵
        % 计算入邻居
        in_neighbors = cell(n, 1);
        for i = 1:n
            in_neighbors{i} = find(obj.A(:, i));
        end
        % 迭代计算SimRank
        for iter = 1:max_iter
            S_old = S;
            for i = 1:n
                for j = 1:n
                    if i == j
                        S(i, j) = 1;
                        continue;
                    end
                    in_i = in_neighbors{i};
                    in_j = in_neighbors{j};
                    if isempty(in_i) || isempty(in_j)
                        S(i, j) = 0;
                        continue;
                    end
                    sum_sim = 0;
                    for a = 1:length(in_i)
                        for b = 1:length(in_j)
                            sum_sim = sum_sim + S_old(in_i(a), in_j(b));
                        end
                    end
                    S(i, j) = C * sum_sim / (length(in_i) * length(in_j));
                end
            end
            if norm(S - S_old, 'fro') < tol
                break;
            end
        end
        scores = S .* obj.train_mask;
    end
    
    function scores = matrix_factorization(obj, dim, lambda, max_iter, learning_rate)
        %MATRIX_FACTORIZATION 矩阵分解方法
        %   dim - 潜在特征维度，默认10
        %   lambda - 正则化参数，默认0.01
        %   max_iter - 最大迭代次数，默认100
        %   learning_rate - 学习率，默认0.01
        if nargin < 2
            dim = 10;
        end
        if nargin < 3
            lambda = 0.01;
        end
        if nargin < 4
            max_iter = 100;
        end
        if nargin < 5
            learning_rate = 0.01;
        end
        [n, ~] = size(obj.A);
        % 初始化用户和物品特征矩阵
        U = randn(n, dim) * 0.01;
        V = randn(n, dim) * 0.01;
        % 获取训练集中的非零元素（即存在的边）
        [rows, cols] = find(obj.train_mask);
        values = ones(length(rows), 1);
        % 随机梯度下降
        for iter = 1:max_iter
            total_error = 0;
            for idx = 1:length(rows)
                i = rows(idx);
                j = cols(idx);
                % 计算预测值和误差
                prediction = U(i, :) * V(j, :)';
                error = values(idx) - prediction;
                total_error = total_error + error^2;
                % 更新特征向量
                U_i_old = U(i, :);
                U(i, :) = U(i, :) + learning_rate * (error * V(j, :) - lambda * U(i, :));
                V(j, :) = V(j, :) + learning_rate * (error * U_i_old - lambda * V(j, :));
            end
            % 添加正则化项
            total_error = total_error + lambda * (norm(U, 'fro')^2 + norm(V, 'fro')^2);
            if mod(iter, 10) == 0
                fprintf('迭代 %d, 误差: %.4f\n', iter, total_error);
            end
        end
        % 计算得分矩阵
        scores = U * V';
        scores = scores .* obj.train_mask;
    end
    
    function [precision, recall, auc] = evaluate(obj, scores, top_k)
        %EVALUATE 评估预测结果
        %   scores - 预测得分矩阵
        %   top_k - 计算precision@k和recall@k的k值
        if nargin < 3
            top_k = 100;
        end
        % 获取测试集中的正样本
        [test_rows, test_cols] = find(obj.test_mask);
        positive_pairs = [test_rows, test_cols];
        num_positives = size(positive_pairs, 1);
        % 获取所有未连接的节点对（负样本+测试正样本）
        negative_mask = (obj.train_mask == 0) & (obj.A == 0) & (eye(size(obj.A)) == 0);
        [negative_rows, negative_cols] = find(negative_mask);
        negative_pairs = [negative_rows, negative_cols];
        % 随机选择与正样本数量相同的负样本
        idx = randperm(size(negative_pairs, 1), num_positives);
        negative_pairs = negative_pairs(idx, :);
        % 合并正负样本
        all_pairs = [positive_pairs; negative_pairs];
        labels = [ones(num_positives, 1); zeros(num_positives, 1)];
        % 获取预测得分
        pred_scores = zeros(size(all_pairs, 1), 1);
        for i = 1:size(all_pairs, 1)
            pred_scores(i) = scores(all_pairs(i, 1), all_pairs(i, 2));
        end
        % 计算AUC
        [~, ~, ~, auc] = perfcurve(labels, pred_scores, 1);
        % 计算Precision@K和Recall@K
        % 获取得分最高的top_k个节点对
        [~, sorted_idx] = sort(pred_scores(1:num_positives), 'descend');
        top_predictions = sorted_idx(1:min(top_k, length(sorted_idx)));
        true_positives = sum(ismember(top_predictions, 1:num_positives));
        precision = true_positives / top_k;
        recall = true_positives / num_positives;
    end
    
    function results = compare_methods(obj, methods, top_k)
        %COMPARE_METHODS 比较不同算法的性能
        %   methods - 要比较的方法列表
        %   top_k - 计算precision@k和recall@k的k值
        if nargin < 2
            methods = obj.methods;
        end
        if nargin < 3
            top_k = 100;
        end
        results = struct();
        for i = 1:length(methods)
            method = methods{i};
            fprintf('正在计算 %s...\n', method);
            try
                % 调用相应的方法
                tic;
                scores = obj.(lower(method))();
                time = toc;
                % 评估性能
                [precision, recall, auc] = obj.evaluate(scores, top_k);
                % 保存结果
                results.(method).scores = scores;
                results.(method).precision = precision;
                results.(method).recall = recall;
                results.(method).auc = auc;
                results.(method).time = time;
                fprintf('%s: Precision@%d=%.4f, Recall@%d=%.4f, AUC=%.4f, 时间=%.2fs\n', ...
                    method, top_k, precision, top_k, recall, auc, time);
            catch ME
                fprintf('计算 %s 时出错: %s\n', method, ME.message);
                results.(method).error = ME.message;
            end
        end
    end
    
    function plot_results(obj, results)
        %PLOT_RESULTS 可视化比较结果
        methods = fieldnames(results);
        num_methods = length(methods);
        precisions = zeros(num_methods, 1);
        recalls = zeros(num_methods, 1);
        aucs = zeros(num_methods, 1);
        times = zeros(num_methods, 1);
        for i = 1:num_methods
            if isfield(results.(methods{i}), 'error')
                continue;
            end
            precisions(i) = results.(methods{i}).precision;
            recalls(i) = results.(methods{i}).recall;
            aucs(i) = results.(methods{i}).auc;
            times(i) = results.(methods{i}).time;
        end
        % 创建图形
        figure('Position', [100, 100, 1200, 800]);
        % 绘制精确度
        subplot(2, 2, 1);
        bar(precisions);
        set(gca, 'XTickLabel', methods, 'XTickLabelRotation', 45);
        title('Precision@K');
        ylabel('Precision');
        grid on;
        % 绘制召回率
        subplot(2, 2, 2);
        bar(recalls);
        set(gca, 'XTickLabel', methods, 'XTickLabelRotation', 45);
        title('Recall@K');
        ylabel('Recall');
        grid on;
        % 绘制AUC
        subplot(2, 2, 3);
        bar(aucs);
        set(gca, 'XTickLabel', methods, 'XTickLabelRotation', 45);
        title('AUC');
        ylabel('AUC');
        grid on;
        % 绘制运行时间
        subplot(2, 2, 4);
        bar(times);
        set(gca, 'XTickLabel', methods, 'XTickLabelRotation', 45);
        title('运行时间');
        ylabel('时间 (秒)');
        grid on;
        % 调整布局
        set(gcf, 'Color', 'w');
    end
end

% 示例使用代码
function example_usage()
    % 生成示例网络（无标度网络）
    n = 100; % 节点数量
    A = create_scale_free_network(n);
    % 创建链路预测对象，使用80%的边作为训练集
    lp = LinkPrediction(A, 0.8);
    % 比较所有方法的性能
    results = lp.compare_methods();
    % 可视化结果
    lp.plot_results(results);
    % 单独使用某个方法
    scores = lp.common_neighbors();
    [precision, recall, auc] = lp.evaluate(scores);
    fprintf('Common Neighbors: Precision=%.4f, Recall=%.4f, AUC=%.4f\n', precision, recall, auc);
end

function A = create_scale_free_network(n)
    %CREATE_SCALE_FREE_NETWORK 生成无标度网络（Barabási-Albert模型）
    %   n - 网络节点数
    % 初始完全图
    m0 = 5; % 初始节点数
    A = zeros(n);
    A(1:m0, 1:m0) = ones(m0) - eye(m0);
    % 添加新节点
    for i = m0+1:n
        % 计算现有节点的度
        degrees = sum(A(1:i-1, 1:i-1), 2);
        total_degree = sum(degrees);
        % 根据度分布选择连接节点
        if total_degree > 0
            prob = degrees / total_degree;
            targets = randsample(1:i-1, m0, true, prob);
        else
            targets = randperm(i-1, min(m0, i-1));
        end
        % 添加连接
        for j = targets
            A(i, j) = 1;
            A(j, i) = 1;
        end
    end
end

% 运行示例
example_usage();

程序说明

这个MATLAB链路预测工具将主流算法集成于同一框架，便于直接对比与调用。具体亮点如下：

1. 核心类 LinkPrediction

一个类统筹所有功能——构造函数读入邻接矩阵，自动按指定比例划分训练集与测试集，随后只需直接调用对应方法即可运行算法。评估与可视化功能也已内置，使用十分便捷。

2. 实现的算法

• Common Neighbors（共同邻居）：两节点共享的邻居越多，未来越可能产生连接。这是最直观的想法，也是许多高级方法的基础。
• Jaccard Coefficient：共同邻居数除以总邻居数，进行归一化处理，有效避免大度节点造成的偏差。
• Adamic-Adar：给予度较小的共同邻居更高权重——“你的朋友不多却愿意与你互动，说明关系更紧密”。
• Preferential Attachment：节点度越大，两节点连接概率越高，符合“富者愈富”的幂律分布特性。
• Katz Index：考虑所有长度的路径，短路径赋予更大权重，通过矩阵求逆一步完成计算。
• Rooted PageRank：从目标节点出发进行随机游走，稳态概率即代表相似度。
• SimRank：若两节点分别连接到相似的节点，则它们自身也相似——一种自洽的迭代定义。
• Matrix Factorization：将网络结构分解为潜在特征，利用SGD进行学习，适用于大规模网络场景。

3. 评估指标

• Precision@K：在预测得分最高的K条边中，实际存在连接的比例。适用于关注“头部推荐”效果的场景。
• Recall@K：所有真实连接中，被Top-K捕捉到的比例，衡量算法对真实连接的覆盖能力。
• AUC：ROC曲线下面积，整体评估排序质量——随机抽取一个正样本与一个负样本，正样本得分高于负样本的概率。

4. 可视化功能

plot_results 方法会生成四个子图：Precision@K、Recall@K、AUC 以及运行时间，结果一目了然。当需要在多种算法之间进行取舍时，这张图能快速提供直观的判断依据。

使用方法

程序末尾已附上完整的示例代码。大致流程为：首先生成一个 Barabási-Albert 无标度网络（100 节点），然后使用 80% 的边作为训练集创建 LinkPrediction 对象，调用 compare_methods 一键运行所有算法并生成图表。也可以单独调用某个方法，例如 common_neighbors，然后自行评估结果。

运行非常简单：将代码复制到 MATLAB 中，直接执行 example_usage 即可。