基于图卷积网络的节点分类

裴来凡

发布于 2024-04-10 14:11:53

830

发布于 2024-04-10 14:11:53

步骤1：下载并加载QM7数据

下载QM7数据集。

dataURL = "http://quantum-machine.org/data/qm7.mat";
outputFolder = fullfile(tempdir,"qm7Data");
dataFile = fullfile(outputFolder,"qm7.mat");

if ~exist(dataFile,"file")
    mkdir(outputFolder);
    disp("Downloading QM7 data...");
    websave(dataFile, dataURL);
    disp("Done.")
end

从MAT文件加载QM7数据。

data = load(dataFile)

data = 包含以下字段的 struct: X: [7165×23×23 single] R: [7165×23×3 single] Z: [7165×23 single] T: [-417.9600 -712.4200 -564.2100 -404.8800 -808.8700 -677.1600 -796.9800 -860.3300 -1.0085e+03 -861.7300 -708.3700 -725.9300 -879.3800 -618.7200 -871.1900 -653.4400 -1.0109e+03 -1.1594e+03 -1.0039e+03 -1.0184e+03 -1.0250e+03 … ] P: [5×1433 int64]

从加载的结构中提取库仑数据和原子序数。

coulombData = double(permute(data.X, [2 3 1]));
atomData = sort(data.Z,2,'descend');

查看第一次观测的原子。

atomData(1,:)

ans = 1×23 single 行向量 6 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

步骤2：预处理图形数据

将训练数据中的库仑矩阵转换为邻接矩阵。

adjacencyData = coulomb2Adjacency(coulombData,atomData);

将图中的前几个分子可视化。

figure
tiledlayout("flow")

for i = 1:9
    % 提取邻接矩阵
    atomicNumbers = nonzeros(atomData(i,:));
    numNodes = numel(atomicNumbers);
    A = adjacencyData(1:numNodes,1:numNodes,i);

    % 将邻接矩阵转换为图
    G = graph(A);

    % 将原子序数转换为符号
    symbols = atomicSymbol(atomicNumbers);

    % 绘图
    nexttile
    plot(G,NodeLabel=symbols,Layout="force")
    title("Molecule " + i)
end

使用直方图可视化每个标签类别的频率。

figure
histogram(categorical(atomicSymbol(atomData)))
xlabel("Node Label")
ylabel("Frequency")
title("Label Counts")

将数据划分为训练、验证和测试分区，分别包含80%、10%和10%的数据。

numObservations = size(adjacencyData,3);
[idxTrain,idxValidation,idxTest] = trainingPartitions(numObservations,[0.8 0.1 0.1]);

adjacencyDataTrain = adjacencyData(:,:,idxTrain);
adjacencyDataValidation = adjacencyData(:,:,idxValidation);
adjacencyDataTest = adjacencyData(:,:,idxTest);

coulombDataTrain = coulombData(:,:,idxTrain);
coulombDataValidation = coulombData(:,:,idxValidation);
coulombDataTest = coulombData(:,:,idxTest);

atomDataTrain = atomData(idxTrain,:);
atomDataValidation = atomData(idxValidation,:);
atomDataTest = atomData(idxTest,:);

预处理训练和验证数据。

[ATrain,XTrain,labelsTrain] = preprocessData(adjacencyDataTrain,coulombDataTrain,atomDataTrain);
size(XTrain)

ans = 88456 1

size(labelsTrain)

ans =

88456 1

[AValidation,XValidation,labelsValidation] = preprocessData(adjacencyDataValidation,coulombDataValidation,atomDataValidation);

使用训练特征的均值和方差对特征进行归一化。

muX = mean(XTrain);
sigsqX = var(XTrain,1);

XTrain = (XTrain - muX)./sqrt(sigsqX);
XValidation = (XValidation - muX)./sqrt(sigsqX);

步骤3：定义深度学习模型

（1）初始化模型参数

创建包含模型参数的结构参数。

parameters = struct;

初始化第一次乘法运算的权重。

numHiddenFeatureMaps = 32;
numInputFeatures = size(XTrain,2);

sz = [numInputFeatures numHiddenFeatureMaps];
numOut = numHiddenFeatureMaps;
numIn = numInputFeatures;
parameters.mult1.Weights = initializeGlorot(sz,numOut,numIn,"double");

初始化第二次乘法运算的权重。

sz = [numHiddenFeatureMaps numHiddenFeatureMaps];
numOut = numHiddenFeatureMaps;
numIn = numHiddenFeatureMaps;
parameters.mult2.Weights = initializeGlorot(sz,numOut,numIn,"double");

初始化第三次乘法运算的权重。

classes = categories(labelsTrain);
numClasses = numel(classes);

sz = [numHiddenFeatureMaps numClasses];
numOut = numClasses;
numIn = numHiddenFeatureMaps;
parameters.mult3.Weights = initializeGlorot(sz,numOut,numIn,"double");

查看参数结构。

parameters

parameters = 包含以下字段的 struct: mult1: [1×1 struct] mult2: [1×1 struct] mult3: [1×1 struct]

查看第一次乘法运算的参数。

parameters.mult1

ans = 包含以下字段的 struct: Weights: [1×32 dlarray]

（2）定义函数模型

定义模型损失函数，指定训练选项。训练1500次，并将Adam优化器的学习率设置为0.01。

numEpochs = 1500;
learnRate = 0.01;

每隔300次验证一次网络。

validationFrequency = 300;

（3）训练模型

初始化Adam的参数。

trailingAvg = [];
trailingAvgSq = [];

将训练和验证功能数据转换为dlarray对象。

XTrain = dlarray(XTrain);
XValidation = dlarray(XValidation);

如果GPU可用，则在GPU上进行训练。

if canUseGPU
XTrain = gpuArray(XTrain);
end

使用onehotencode函数将训练和验证标签转换为一个热编码向量。

TTrain = onehotencode(labelsTrain,2,ClassNames=classes);
TValidation = onehotencode(labelsValidation,2,ClassNames=classes);

初始化TrainingProgressMonitor对象。

monitor = trainingProgressMonitor( ...
    Metrics=["TrainingLoss","ValidationLoss"], ...
    Info="Epoch", ...
    XLabel="Epoch");

groupSubPlot(monitor,"Loss",["TrainingLoss","ValidationLoss"])

使用自定义训练循环训练模型。

epoch = 0;

while epoch < numEpochs && ~monitor.Stop
epoch = epoch + 1;

    % 评估模型损失函数和梯度下降
    [loss,gradients] = dlfeval(@modelLoss,parameters,XTrain,ATrain,TTrain);

    % 使用Adam优化器更新网络参数
    [parameters,trailingAvg,trailingAvgSq] = adamupdate(parameters,gradients, ...
trailingAvg,trailingAvgSq,epoch,learnRate);

    % 训练损失和迭代次数
recordMetrics(monitor,epoch,TrainingLoss=loss);
updateInfo(monitor,Epoch=(epoch+" of "+numEpochs));

    % 验证指标
if epoch == 1 || mod(epoch,validationFrequency) == 0
YValidation = model(parameters,XValidation,AValidation);
lossValidation = crossentropy(YValidation,TValidation,DataFormat="BC");

    % 验证损失
recordMetrics(monitor,epoch,ValidationLoss=lossValidation);
end

monitor.Progress = 100*(epoch/numEpochs);
end

（4）测试模型

使用测试数据测试模型。

[ATest,XTest,labelsTest] = preprocessData(adjacencyDataTest,coulombDataTest,atomDataTest);
XTest = (XTest - muX)./sqrt(sigsqX);

将测试功能数据转换为dlarray对象。

XTest = dlarray(XTest);

对数据进行预测，并使用onehotdecode函数将概率转换为分类标签。

YTest = model(parameters,XTest,ATest);
YTest = onehotdecode(YTest,classes,2);

计算精度。

accuracy = mean(YTest == labelsTest)

accuracy = 0.9001

使用混淆图函数计算混淆矩阵。

figure
cm = confusionchart(labelsTest,YTest, ...
    ColumnSummary="column-normalized", ...
    RowSummary="row-normalized");
title("GCN QM7 Confusion Chart");

（5）使用新数据进行预测

对未标记的数据进行一些预测。

numObservationsNew = 4;
adjacencyDataNew = adjacencyDataTest(:,:,1:numObservationsNew);
coulombDataNew = coulombDataTest(:,:,1:numObservationsNew);

predictions = modelPredictions(parameters,coulombDataNew,adjacencyDataNew,muX,sigsqX,classes);

在绘图中可视化预测。

figure
tiledlayout("flow")

for i = 1:numObservationsNew
    % 提取未添加的数据
    numNodes = find(any(adjacencyDataTest(:,:,i)),1,"last");

    A = adjacencyDataTest(1:numNodes,1:numNodes,i);

    % 绘图
    nexttile
    G = graph(A);
    plot(G,NodeLabel=string(predictions{i}),Layout="force")
    title("Observation " + i + " Prediction")
end

（6）支持功能

预处理数据功能。

function [adjacency,features,labels] = preprocessData(adjacencyData,coulombData,atomData)

[adjacency, features] = preprocessPredictors(adjacencyData,coulombData);
labels = [];

% 将标签转换为类别标签
for i = 1:size(adjacencyData,3)
    % 提取并添加未添加的数据
    T = nonzeros(atomData(i,:));
    labels = [labels; T];
end

labels2 = nonzeros(atomData);
assert(isequal(labels2,labels2))

atomicNumbers = unique(labels);
atomNames =  atomicSymbol(atomicNumbers);
labels = categorical(labels, atomicNumbers, atomNames);

end

预处理预测器函数。

function [adjacency,features] = preprocessPredictors(adjacencyData,coulombData)

adjacency = sparse([]);
features = [];

for i = 1:size(adjacencyData, 3)
% 提取未添加的数据
numNodes = find(any(adjacencyData(:,:,i)),1,"last");

A = adjacencyData(1:numNodes,1:numNodes,i);
X = coulombData(1:numNodes,1:numNodes,i);

% 从库仑矩阵的对角线中提取特征向量
X = diag(X);

% 添加提取的数据
adjacency = blkdiag(adjacency,A);
features = [features; X];
end

end

函数模型。

function Y = model(parameters,X,A)

ANorm = normalizeAdjacency(A);

Z1 = X;

Z2 = ANorm * Z1 * parameters.mult1.Weights;
Z2 = relu(Z2) + Z1;

Z3 = ANorm * Z2 * parameters.mult2.Weights;
Z3 = relu(Z3) + Z2;

Z4 = ANorm * Z3 * parameters.mult3.Weights;
Y = softmax(Z4,DataFormat="BC");

end

模型损失函数。

function [loss,gradients] = modelLoss(parameters,X,A,T)

Y = model(parameters,X,A);
loss = crossentropy(Y,T,DataFormat="BC");
gradients = dlgradient(loss, parameters);

end

模型预测函数。

function predictions = modelPredictions(parameters,coulombData,adjacencyData,mu,sigsq,classes)

predictions = {};
numObservations = size(coulombData,3);

for i = 1:numObservations
% 提取未添加的数据
numNodes = find(any(adjacencyData(:,:,i)),1,"last");
A = adjacencyData(1:numNodes,1:numNodes,i);
X = coulombData(1:numNodes,1:numNodes,i);

% 预处理数据
[A,X] = preprocessPredictors(A,X);
X = (X - mu)./sqrt(sigsq);
X = dlarray(X);

% 预测
Y = model(parameters,X,A);
Y = onehotdecode(Y,classes,2);
predictions{end+1} = Y;
end

end

归一化邻接矩阵函数。

function ANorm = normalizeAdjacency(A)

% 将自连接添加到邻接矩阵
A = A + speye(size(A));

% 计算度数的平方倒数
degree = sum(A, 2);
degreeInvSqrt = sparse(sqrt(1./degree));

% 归一化邻接矩阵
ANorm = diag(degreeInvSqrt) * A * diag(degreeInvSqrt);

end

链接：http://quantum-machine.org/data/qm7.mat

https://ww2.mathworks.cn/help/deeplearning/ug/node-classification-using-graph-convolutional-network.html?s_tid=srchtitle_site_search_10_LORENZ

本文参与腾讯云自媒体分享计划，分享自微信公众号。

原始发表：2024-04-01，如有侵权请联系 cloudcommunity@tencent.com 删除

测试