图像识别3：梯度下降和LDA线性模型实验

图像识别3：梯度下降和LDA线性模型实验

• 写在最前面
• 一、实验内容
• 二、实验结果
• 三、实验源码

• 实验六：梯度下降一元回归
• 实验七：LDA一元回归

一、实验内容

1. 根据梯度下降法、闭式解完成一元线性回归实验，并比较两种解下的实验结果。
   
   

划分训练集和测试集比例(1:2 划分）。通过 randperm 函数得到一个打乱序列的 b，该序列包含（1，n），本文中 n 为合并后矩阵的行数。由于题目要求的 1:2 划分容易出现不完美分割的情况，因此对划分的数量进行取整操作。然后 c,d 两个序列分别取打乱序列的前 33%、后 66%，根据 c,d 两个序列的数值取原矩阵对应的列数，最后将按顺序将每次划分结果对应保存于对应的训练集和测试集中，并分别绘制训练集和测试集的散点图。

遍历测试集与训练集，计算两者两两之间的欧式距离。绘图展示梯度下降线性回归，最小二乘线性回归的图像。

4. 根据训练数据，学习投影矩阵。然后将 LDA 在训练样本上的低维表示结果可视化。最后使用距离最短对测试样本进行分类。
   首先定义线性判别回归函数 LDA。其中，mappedX 为低维数据降维后的坐标，mapping 为映射的信息，M 为，X 为样本矩阵，label 为样本标签，no_dims: 投影空间(低维空间)的维度。

初始化。①设置特征默认维数为 2，提高函数的复用性。然后通过②bsxfun 函数对函数进行标准化，优点为当 A 和 B 的维度不相等，并且 A 和 B各自有一个维度为 1 时，函数自动调整维度，目的也是提高函数复用性。之后③去掉矩阵标签中重复的元素，将标签类别保存于 unique 中。最后初始化用到的矩阵。

通过①计算总协方差矩阵（全局散度矩阵）St = np.dot((X - mean).T, X - mean)，②St = np.dot((X - mean).T, X - mean)，③全局散度矩阵 St=Sw+Sb，所以 Sb = St-Sw 三个式子计算类内离散度矩阵 Sw。同时确保降维后的特征点维度小于 max，提高代码的安全性。

按降序排列特征值和特征向量。M 变换矩阵由 v 的最大的 nc-1 个特征值所对应的特征向量构成。对 inv(Sw)*Sb 进行特征分解[V,D]=eig(A,B) eig(A,B)返回方阵 A 和 B 的 N 个广义特征值，构成 N×N 阶对角阵 D， 其对角线上的N 个元素即为相应的广义特征值，同时将返回相应的特征向量构成 N×N 阶满秩矩阵，且满足 AV=BVD。之后判断数组的元素是否是 NaN，并提取前no_dims 个特征向量。

最后计算投影后的（二维）坐标值。

主函数部分，先通过 xlsread 或 csread 函数读取 csv 文件数据。

在 LDA 多分类后，将返回的四维降维成的二维坐标值反应到图像上，通过不同颜色和形状区分不同的标签，计算坐标两两之间的欧氏距离。

二、实验结果

将其改为 length 函数。

梯度下降得到的 theta 的值:

w:1.159369,

b:-3.369070

梯度下降法一元回归的欧式距离为 591.925823,

闭式解为 594.443011

可以得到结果，梯度下降和最小二乘线性回归的直线几乎重合，差别不大。
训练集[日本，中国]的标签分别为 1,2。
结果很合适，可以看到点就在散点图聚集处。

三、实验源码

实验六：梯度下降一元回归

% 1.根据梯度下降法完成一元线性回归实验。
% 2.根据闭式解完成一元线性回归实验。
% 3.比较两种解下的实验结果。

% 定义多元线性回归计算代价损失函数
% 参数进行线性回归拟合X和y中的数据点
function J = computeCost(X, y, theta) %计算使用theta作为的成本
	m = length(y); % number of training examples
	J = 0; % 返回正确的变量
    % J <= 10e-5; %保证误差小于10的-5次方
	predictions = X * theta;
	J = 1/(2*m)*(predictions - y)'*(predictions - y);
    % J = 1/(2*m)*(predictions - y).^2; %线性函数的求解公式
end


%定义一个实现梯度下降的函数
function [theta, J_history] = gradientDescent(X, y, theta, alpha, num_iters)
m = length(y);%取数据的长度
J_history = zeros(num_iters, 1);%定义J_history为num_iters行1列的向量
%其中，num_iters是在调用该梯度下降函数时的参数
for iter = 1:num_iters
    S = (1 / m) * (X' * (X * theta - y));%相当于求导
    theta = theta - alpha .* S;%theta的更新
    J_history(iter) = computeCost(X, y, theta);%在这里调用computeCost(X, y, theta)
 %相当于在慢慢的减小代价函数
end
end


% 主函数
% 1.根据梯度下降法完成一元线性回归实验。
% 2.根据闭式解完成一元线性回归实验。
% 3.比较两种解下的实验结果。
clc;clear;
close all;
%% 数据导入
data = load('D:\Desktop\大三上\神经网络\实验6\ex1data1.txt'); %读入文件数据
x=data(:,1); %取data中第一列，人口每行一个样本
y=data(:,2); %取data中第二列，利润每行一个样本
% plotData(x,y); %显示样本数据图
%% 初始化值
m=length(x); %取x的长度，即数据的个数
x=[ones(m,1),data(:,1)]; %将原来m*1的x，多拓展了一个元素列
%全为1的一列，变为m*2的矩阵，其中的一列为全1列，主要是为了方便后面的运算
theta=zeros(2,1); %给theta先赋值为0，此theta相当于theta0和theta1
%组成的2*1的向量
iterations=1500; %设迭代次数为1500
alpha=0.01; %设梯度下降步长为0.01
%% 划分训练集和测试集比例(1:2划分）
b=randperm(length(x(:,1))); %打乱列序列
e=round(length(x(:,1))*1/3);
c=b(1:e); %取打乱序列的前33%
d=b(e+1:end); %取打乱序列的后66% %end直到数组结尾

atrainx=x(c,:); %(:,:)为取原矩阵行和列，(:,1:50)为取原矩阵列，行取前33%
atrainy=y(c,:);
atestx=x(d,:); %(:,d)为取原矩阵列，行取剩余的随机66%的行
atesty=y(d,:);

plot(atrainx(:,2),atrainy,'ro','MarkerSize',10) %显示了训练数据样本数据散点图
hold on; %图保留
plot(atestx(:,2),atesty,'yo','MarkerSize',10) %显示了测试数据样本数据散点图
hold on; %图保留
%% 梯度下降
theta=gradientDescent(atrainx,atrainy,theta,alpha,iterations); %调用梯度
%下降更新theta
% matlab索引从1开始
fprintf('梯度下降得到的theta的值:\nw:%f,\nb:%f\n',theta(2),theta(1)) %theta(2)为b
%% 闭式解
p = polyfit(atrainx(:,2),atrainy,1);
%% 计算测试集与训练集两者之间的欧式距离
predicty=zeros(2,length(atesty)); %初始化两列预测y的数据集
% for i=1:length(d)
for i=1:d
    predicty(1,i)=[1,i]*theta; %第一列预测y为梯度下降算法求得的解
    predicty(2,i)=polyval(p,i); %第一列预测y为闭式解求得的解
    % predicty(2,i)=i*p(1)+p(2); %polyval(p,i)就是i*p(1)+p(2)的意思
end
distance1=sqrt(sum(abs(predicty(1,:)-atesty)).^2);  %sum(x.^2)计算数组的平方和
distance2=sqrt(sum(abs(predicty(2,:)-atesty)).^2);  %sum(x.^2)计算数组的平方和
fprintf('梯度下降法一元回归的欧式距离为 %f,\n闭式解为 %f\n',mean(distance1),mean(distance2));
%% 绘图展示结果
hold on; %保持图像不被刷新
plot(atrainx(:,2),atrainx*theta,'b-','MarkerSize',10) %显示了梯度下降拟合的直线
hold on; %保持图像不被刷新
plot(atrainx(:,2),atrainx(:,2)*p(1)+p(2),'b-','MarkerSize',10) %显示了最小二乘拟合的直线
legend('训练数据','测试数据','梯度下降线性回归','最小二乘线性回归') %图例标注

实验七：LDA一元回归

clear;clc
%训练数据集数据和标签
%% 读取路径下表一数据
%方法一：xlsread
% M = xlsread('D:\Desktop\实验七\EduRecivedFiles\LDA数据.csv',1,'C2:G23');
%方法二：csvread
% 文件名，开始读取位置的行号（row）和列号（col）注意从0开始
% range = [R1 C1 R2 C2]，这里（R1，C1）是读取区域的左上角，（R2，C2）是读取区域的右下角
M1 = csvread('D:\Desktop\实验七\EduRecivedFiles\LDA数据.csv', 1, 2, [1 2 22 6]);
x=M1(:,1:end-1);
y=M1(:,end);
%% LDA多分类
[mappedX,mapping,W]=FisherLDA2(x,y,2);
%% 结果可视化
figure
temp1=find(y==1);
plot(mappedX(temp1,1),mappedX(temp1,2),'r^','MarkerSize',18)
hold on
temp2=find(y==2);
plot(mappedX(temp2,1),mappedX(temp2,2),'bo','MarkerSize',18)
hold on
temp3=find(y==3);
plot(mappedX(temp3,1),mappedX(temp3,2),'black*','MarkerSize',18)
hold on
% Z=-4:5;
% plot(Z,Z*W,'-k')
% axis([-4 5 -4 2]);
%% 测试样本分类
% 绘制测试样本
temp0=find(y==0);
plot(mappedX(21,1),mappedX(21,2),'y.','MarkerSize',18)
plot(mappedX(22,1),mappedX(22,2),'green.','MarkerSize',18)
legend('标签1','标签2','标签3','测试样本1日本','测试样本2中国');
hold on
% 计算两者之间的欧式距离
tempdis1=inf;tempdis2=inf;  %初始化最短距离
tempj=0;%初始化最短距离点的坐标
for i=1:length(mappedX)-2
    distance1=sqrt(sum(abs(mappedX(i,:)-mappedX(21,:))).^2);  %sum(x.^2)计算数组的平方和
    distance2=sqrt(sum(abs(mappedX(i,:)-mappedX(22,:))).^2);  %sum(x.^2)计算数组的平方和
    if(distance1<tempdis1)
        tempdis1=distance1;tempj1=i; % 如果距离最小则判断为同一类型
        res1=y(i,:); % 将最短距离的训练集标签赋值给测试集
    end
    if(distance2<tempdis2)
        tempdis2=distance2;tempj2=i; % 如果距离最小则判断为同一类型
        res2=y(i,:); % 将最短距离的训练集标签赋值给测试集
    end
end
fprintf('训练集[日本，中国]的标签分别为%d,%d\n',res1,res2)
% 绘制一维图像()
% [mappedX1,mapping1,W1]=FisherLDA2(x,y,1);
% Z=-4:5;
% plot(Z,Z*W(1),'-k')
% axis([-4 5 -4 2]);
% hold on
% legend('标签1','标签2','标签3','测试样本','一维直线');

function [mappedX, mapping, M] = FisherLDA2(X, labels, no_dims)
% X:        n*d matrix,each row is a sample;
% labels:   n*k matrix,each is the class label
% no_dims:  投影空间(低维空间)的维度(no_dimms >= 1, default = 2，max = len(labels)-1
% mappedX： 低维数据降维后的坐标
% mapping： 映射的信息

%% 初始化
% 设置特征点默认的维数为2
if ~exist('no_dims', 'var') || isempty(no_dims)
    no_dims = 2;
end

% 标准化数据集
mapping.mean = mean(X, 1);
% 当A和B的维度不相等，并且A和B各自有一个维度为1时，函数自动调整维度
X = bsxfun(@minus, X, mapping.mean); % 减法操作

% 去掉矩阵标签中重复的元素
[classes, bar, labels] = unique(labels);
nc = length(classes);

% 初始化类内离散度矩阵 Sw
Sw = zeros(size(X, 2), size(X, 2));
% 计算总协方差矩阵（全局散度矩阵）St = np.dot((X - mean).T, X - mean)
St = cov(X);

%% 全局散度矩阵St=Sw+Sb，so Sb = St-Sw
% St = np.dot((X - mean).T, X - mean)
% 计算类内离散度矩阵Sw
for i=1:nc
    cur_X = X(labels == i,:);
    % 更新类内离散度矩阵Sw
    C = cov(cur_X);
    p = size(cur_X, 1) / (length(labels) - 1);
    Sw = Sw + (p * C);
end

% 计算类间离散度矩阵Sb = St-Sw
Sb = St - Sw;
Sb(isnan(Sb)) = 0; Sw(isnan(Sw)) = 0;
Sb(isinf(Sb)) = 0; Sw(isinf(Sw)) = 0;

% 确保嵌入特征空间中特征点的维数小于max
if nc <= no_dims
    no_dims = nc - 1;
    warning(['目标维度降为： ' num2str(no_dims) ]);
end

%% M变换矩阵由v的最大的nc-1个特征值所对应的特征向量构成
% 对inv(Sw)*Sb进行特征分解[V,D]=eig(A,B)
% eig(A,B)返回方阵A和B的N个广义特征值，构成N×N阶对角阵D，
% 其对角线上的N个元素即为相应的广义特征值，同时将返回相应的特征向量构成N×N阶满秩矩阵，且满足AV=BVD。
[M, lambda] = eig(Sb, Sw);

% 按降序排列特征值和特征向量
lambda(isnan(lambda)) = 0; % 判断数组的元素是否是NaN
[lambda, ind] = sort(diag(lambda), 'descend'); %diag(lambda)提取lambda矩阵的对角线
M = M(:,ind(1:min([no_dims size(M, 2)]))); %提取前no_dims个特征向量，ind——最小值下标

%% 计算投影后的中心值
mappedX = X * M;

% 存储测试集投影中心值
mapping.M = M;
mapping.val = lambda;