存档在 ‘模式识别&机器学习’ 分类

矩阵的迹 特征值

2017年11月25日

矩阵的迹(trace)

X∈P(n×n,X=(xii)的主对角线上的所有元素之和称之为X的迹,记为tr(X),即tr(X)=∑xii

性质:

(1)

设有N阶矩阵A,那么矩阵A的迹(用tr(A)表示)就等于A的特征值的总和,也即A矩阵的主对角线元素的总和。

1.迹是所有对角元的和

2.迹是所有特征值的和

3.某些时候也利用tr(AB)=tr(BA)来求迹

(2)

奇异值分解(Singular value decomposition )

奇异值分解非常有用,对于矩阵A(p*q),存在U(p*p),V(q*q),B(p*q)(由对角阵与增广行或列组成),满足A = U*B*V

U和V中分别是A的奇异向量,而B是A的奇异值。AA’的特征向量组成U,特征值组成B’B,A’A的特征向量组成V,特征值(与AA’相同)组成BB’。因此,奇异值分解和特征值问题紧密联系。

如果A是复矩阵,B中的奇异值仍然是实数。

SVD提供了一些关于A的信息,例如非零奇异值的数目(B的阶数)和A的阶数相同,一旦阶数确定,那么U的前k列构成了A的列向量空间的正交基。

 

 

矩阵的特征值(eigenvalue)

设 A 是n阶方阵,如果存在数m和非零n维列向量 x,使得 Ax=mx 成立,则称 m 是A的一个特征值(characteristic value)或本征值(eigenvalue)。

求解矩阵特征值的方法:

Ax=mx,等价于求m,使得(mE-A)x=0,其中E是单位矩阵,0为零矩阵。

|mE-A|=0,求得的m值即为A的特征值。|mE-A| 是一个n次多项式,它的全部根就是n阶方阵A的全部特征值,这些根有可能相重复,也有可能是复数。

如果n阶矩阵A的全部特征值为m1 m2 … mn,则|A|=m1*m2*…*mn

同时矩阵A的迹是特征值之和:tr(A)=m1+m2+m3+…+mn[1]

如果n阶矩阵A满足矩阵多项式方程g(A)=0, 则矩阵A的特征值m一定满足条件g(m)=0;特征值m可以通过解方程g(m)=0求得

机器学习: Canonical Correlation Analysis 典型相关分析

2017年10月16日

Canonical Correlation Analysis(CCA)典型相关分析也是一种常用的降维算法。我们知道,PCA(Principal Component Analysis) 主分量分析将数据从高维映射到低维空间同时,保证了数据的分散性尽可能地大, 也就是数据的方差或者协方差尽可能大。而LDA(Linear Discriminant Analysis) 线性判别分析则利用了类标签,利用一种监督学习的方法,将数据从高维空间映射到低维空间时,让不同类的数据尽可能地分开而同一类的数据尽可能地聚合。

但是,有的时候,想知道多个线性空间之间的相关性。比如有的时候我们会从图像中提取各种特征,每一种特征都可以构成一个线性空间,为了分析这些空间之间的相关性,我们可以利用CCA 来做分析。
假设我们有两个特征空间,\(S1={\mathbf{x}_{1} \in R^{d1}}\)\(S2={\mathbf{x}_{2} \in R^{d2}}\)

\(\mathbf{x} = \begin{pmatrix}
\mathbf{x}_{1} \\
\mathbf{x}_{2}
\end{pmatrix}
\quad E(\mathbf{x}) = \begin{pmatrix}
\mathbf{\mu}_{1} \\
\mathbf{\mu}_{2}
\end{pmatrix}
\quad \Sigma = \begin{pmatrix}
\Sigma_{11} & \Sigma_{12} \\
\Sigma_{21} & \Sigma_{22}
\end{pmatrix}\)

可以看到,\(\Sigma_{12}=\Sigma_{21}^{T}\),\(\Sigma\) 称为协方差矩阵。我们引入投影向量 \(\mathbf{a}\), \(\mathbf{b}\), 假设投影之后的变量满足:

\(u=\mathbf{a}^{T} \mathbf{x}_{1} \quad v=\mathbf{b}^{T} \mathbf{x}_{2}\)

可以进一步算出 \(u,v\)的方差和协方差:

\(\text{var}(u)= \mathbf{a}^{T} \Sigma_{11} \mathbf{a}, \quad \text{var}(v)=\mathbf{b}^{T} \Sigma_{21} \mathbf{b}, \quad cov(u,v)=\mathbf{a}^{T} \Sigma_{12} \mathbf{b}\)

可以计算出 \(u,v\)的相关系数:

\(Corr(u,v)=\frac{\text{cov}(u,v)}{\sqrt{\text{var}(u)} \sqrt{\text{var}(v)}}\)

将u,v的表达式代入,可以得到:

\(Corr(u,v)=\frac{\mathbf{a}^{T} \Sigma_{12} \mathbf{b}}{\sqrt{\mathbf{a}^{T} \Sigma_{11} \mathbf{a}} \sqrt{\mathbf{b}^{T} \Sigma_{22} \mathbf{b}}}\)

我们的目标是让相关系数\(Corr(u,v)\)尽可能地大。为了求解\(\mathbf{a}\), \(\mathbf{b}\), 可以固定分母而让分子最大化,所以上面的函数可以变成:

\(\max_{\mathbf{a}, \mathbf{b}} \mathbf{a}^{T} \Sigma_{12} \mathbf{b}\)

构造拉格朗日等式:

\(L=\mathbf{a}^{T} \Sigma_{12} \mathbf{b}-\frac{\lambda_{1}}{2}(\mathbf{a}^{T} \Sigma_{11} \mathbf{a}-1)-\frac{\lambda_{2}}{2}(\mathbf{b}^{T} \Sigma_{22} \mathbf{b}-1)\)

L 分别对\(\mathbf{a}\), \(\mathbf{b}\)求导,可以得到:

\(\frac{\partial L}{\partial \mathbf{a}}= \Sigma_{12} \mathbf{b}- \lambda_{1} \Sigma_{11} \mathbf{a}=0\)

\(\frac{\partial L}{\partial \mathbf{b}}= \Sigma_{21} \mathbf{a}- \lambda_{2} \Sigma_{22} \mathbf{b}=0\)

根据约束条件,可以得到:

\(\lambda_{1}=\lambda_{2}=\mathbf{a}^{T} \Sigma_{12} \mathbf{b}\)

所以只要求出\(\lambda_{1}\) 或者\(\lambda_{2}\)就可以得到最大的相关系数。令\( \lambda=\lambda_{1}=\lambda_{2}\).
通过上面的偏导数,我们可以得到:

\(\Sigma_{11}^{-1}\Sigma_{12} \mathbf{b}= \lambda \mathbf{a}\)

\(\Sigma_{22}^{-1}\Sigma_{21} \mathbf{a}= \lambda \mathbf{b}\)

写成矩阵形式:

\(\begin{pmatrix}
\Sigma_{11}^{-1} & 0 \\
0 & \Sigma_{22}^{-1}
\end{pmatrix} \begin{pmatrix}
0 & \Sigma_{12} \\
\Sigma_{21} & 0
\end{pmatrix} \begin{pmatrix}
\mathbf{a} \\
\mathbf{b}
\end{pmatrix}=\lambda \begin{pmatrix}
\mathbf{a} \\
\mathbf{b}
\end{pmatrix}\)

令:

\(B= \begin{pmatrix}
\Sigma_{11} & 0 \\
0 & \Sigma_{22} \end{pmatrix}, \quad A= \begin{pmatrix}
0 & \Sigma_{12} \\
\Sigma_{21} & 0
\end{pmatrix} \quad \mathbf{w}=\begin{pmatrix}
\mathbf{a} \\
\mathbf{b}
\end{pmatrix}\)

那么,上式可以表示成:

\(B^{-1}A\mathbf{w}=\lambda \mathbf{w}\)

所以,\(\lambda \)和 \(\mathbf{w}\) 就是\(B^{-1}A\) 的特征值和特征向量。我们可以求出 \(B^{-1}A\)的特征值和特征向量,然后利用特征向量将原来的特征
\(\mathbf{x}_{1}, \mathbf{x}_{2}\)做映射。对应特征值 \(\lambda \)的求解,可以有更简单的方法,从上面的偏导数,我们可以得到如下等式:

\(\Sigma_{11}^{-1}\Sigma_{12}\Sigma_{22}^{-1}\Sigma_{21} \mathbf{a}=\lambda^{2} \mathbf{a}\)

我们可以利用上面的表达式求出\(\lambda \) 和 \(\mathbf{a}\),然后再待会上面的偏导数等式求出 \(\mathbf{b}\).
\(\lambda \)就是 \(u,v\)的相关系数,\(u,v\) 就是一对典型变量(canonical variables)。按照 \(B^{-1}A\) 的特征值从大到小排列,可以求出一系列的典型变量。特征值越大,说明典型变量的相关性越强。
Reference:
http://blog.csdn.net/wenyusuran/article/details/35290683

PCA and Whitening

2017年10月14日

文献:http://ufldl.stanford.edu/wiki/index.php/Exercise:PCA_and_Whitening

(一):数据准备

UFLDL下载的文件中,包含数据集IMAGES_RAW,它是一个512*512*10的矩阵,也就是10幅512*512的图像,可以看下其中的第二幅,使用extract_show函数。

function extract_show(IMAGESr, index)
%load IMAGES_RAW;
IMAGES_ = IMAGESr;
imshow(IMAGES_(:,:,index));

FE957034-C785-4284-B361-2B0100D32847

(1)载入数据

利用sampleIMAGESRAW函数,从IMAGES_RAW中提取numPatches个图像块,每个图像块大小为patchSize,并将提取到的图像块按列存放,分别存放在在矩阵patches的每一列中,即patches(:,i)存放的是第i个图像块儿的所有像素值。如下为产生的原始图像部分。
1

(2)数据去均值化处理

将每一个图像块儿的所有像素值都减去该图像块儿的平均像素值,实现数据的去均值化。
这里要注意一下第三行。理论上应该使用第二行做均值归0,但是图像每个像素点之间是很相关的。

[dim,sampleN] = size(x);%维数与样本数
%x = x - repmat(mean(x,2),1,sampleN);%标准PCA均值归0
x = x - repmat(mean(x,1),dim,1);

2

(二)PCA处理

该部分分为两部分:

(1)进行PCA计算

这里仅仅对数据x进行旋转得到xrot,而不进行主成分的提取,具体地:
①计算数据x的协方差矩阵sigma
②对sigma进行特征分解,利用matlab的eig函数,从而得到sigma的特征向量构成的矩阵U

[U,S,V]=eig(sigma);

U=[u1,…,ui,…,un],它的每一列分别是sigma的特征向量,n是输入数据的特征维数,S=diag([λ1,…λi,…,λn])是由sigma的特征值作为对角元素的对角阵,ui和λi相对应,为了后续的计算,这里要将U的各列次序进行调换,使得调换后的各列所对应的特征值大小依次递减,调换后的矩阵仍记作U,相应的特征值对角阵仍即为S,即:

U=[u1,...,ui,...,un],S=diag([λ1,...λi,...,λn])

满足:

λ1>=...>=λi>=...>=λn

③利用矩阵U对数据x进行旋转,得到xrot,即

xrot=U‘*x

完整代码如下:

xRot = zeros(size(x)); % You need to compute this 144 * 10000
S = x*x'/sampleN; %协方差矩阵
[U,D] = eig(S); %特征向量与特征值
U = fliplr(U);%U为列特征向量,D为特征值对角阵
D = fliplr(fliplr(D)'); %按照特征值由大到小排列
xRot = U'*x;

旋转变换后图为
3

(2)对旋转后的数据求解协方差矩阵covar

将协方差矩阵可视化,观察得到的选择后的数据是否正确。PCA保证选择后的数据的协方差矩阵是一个对角阵,如果covar是正确的,那么它的图像应该是一个蓝色背景,并且在对角线位置有一斜线。这里显示协方差矩阵covar利用了matlab的imagesc,imagesc(covar)的作用是:把矩阵covar以图像形式显示出来,矩阵中不同的数值会被赋予不同的颜色,得到的协方差矩阵的图像如下:可以看到,图像除了对角线位置外,其余部分颜色都相同。

covar = zeros(size(x, 1)); % You need to compute this
% Visualise the covariance matrix. You should see a line across the
% diagonal against a blue background.
covar = xRot*xRot';%获得变换后数据的协方差矩阵

4

(三)获取满足条件的主成分

本部分,找到满足条件的主成分的个数k,也就是找到最小的k值,使得(λ1+…+ λk)/(λ1+…+ λn)>某个百分数,如99%

k = 0; % Set k accordingly
var_total = sum(D(:))               %总方差
var_cur = 0;                        %方差计数
for k = 1:dim
    var_cur = var_cur + D(k,k);     %累计
    if (var_cur / var_total) > 0.90 %如果能量已>90%,则跳出for语句
        break;
    end
end

(四)利用找到的主成分个数,对数据进行降维

在第三步,已经找到了数字k,也就是,保留数据的k个主成分就满足了要求。在该步,将对数据x进行降维,只留下k个主成分。同时,为了观察降维后的数据的好坏,在利用U(:,k)将降维后的数据变换会原来的维数,也就是得到了原数据的近似恢复数据。并利用网格将恢复出的图像显示出,与原图像进行比较,下面第一幅图是由降维后的图像恢复出的原数据,下图是相应的原数据,可以发现,降维后的数据基本可以恢复出于原数据非常相近的数据。

xHat = zeros(size(x));  % You need to compute this
U_reduce = zeros(dim);          %初始化
U_reduce(:, 1:k) = U(:, 1:k);   %去掉特征值小的特征向量(降维)
y = U_reduce'*x;                %变换
xHat = U_reduce*y;              %反变换回来

% Visualise the data, and compare it to the raw data
% You should observe that the raw and processed data are of comparable quality.
% For comparison, you may wish to generate a PCA reduced image which
% retains only 90% of the variance.

figure('name',['#5 PCA processed images ',sprintf('(%d / %d dimensions)', k, size(x, 1)),'']);
display_network(xHat(:,randsel));
figure('name','Raw images');
display_network(x(:,randsel));

5
1
假如保存前90%的能量,PCA可以用44个特征代表144个特征。

(五)PCA白化+正则化

该部分分为两步

(1)执行具有白化和正则化的PCA

首先,对数据进行旋转(利用特征矩阵U),然后,利用特征值对旋转后的数据进行缩放,实现白化。同时,在利用特征值缩放时,利用参数ε对特征值进行微调,实现正则化。

epsilon = 0.1;
xPCAWhite = zeros(size(x));

eig_value = diag(D);                            %特征值
normalize = sqrt(eig_value+epsilon);            %白化系数
xPCAWhite = U'*x ./ (normalize*ones(1,sampleN));%白化

271631124694609

(2)计算白化后的数据的协方差矩阵,观察该协方差矩阵

如果加入了正则化项,则该协方差矩阵的对角线元素都小于1,如果没有加入正则项(即仅有旋转+白化),则该协方差矩阵的对角线元素都为1(实际上,是令ε为一个极小的数),下图是白化后数据的协方差矩阵对应的图像,上图是加入正则化后的结果,下图是没有加入正则化后的结果。

covar = xPCAWhite*xPCAWhite' / sampleN; %协方差矩阵 
figure('name','Visualisation of covariance matrix');
imagesc(covar);

8

(六)ZCA白化

ZCA白化,就是在PCA白化的基础上做了一个旋转,即

xZCAWhite = zeros(size(x));
xZCAWhite = U * xPCAWhite;  %ZCA白化 

271923492655269
下面的第一幅图是ZCA白化后的结果图,第二幅图是相应的原始图像。可以看到,ZCA白化的结果图似乎是原始图像的边缘。
white
1

Stanford UFLDL教程 MATLAB Modules

2017年10月14日

MATLAB Modules

Sparse autoencoder |sparseae_exercise.zip

  • checkNumericalGradient.m – Makes sure that computeNumericalGradient is implmented correctly
  • computeNumericalGradient.m – Computes numerical gradient of a function (to be filled in)
  • display_network.m – Visualizes images or filters for autoencoders as a grid
  • initializeParameters.m – Initializes parameters for sparse autoencoder randomly
  • sampleIMAGES.m – Samples 8×8 patches from an image matrix (to be filled in)
  • sparseAutoencoderCost.m -Calculates cost and gradient of cost function of sparse autoencoder
  • train.m – Framework for training and testing sparse autoencoder

Using the MNIST Dataset |mnistHelper.zip

  • loadMNISTImages.m – Returns a matrix containing raw MNIST images
  • loadMNISTLabels.m – Returns a matrix containing MNIST labels

PCA and Whitening |pca_exercise.zip

  • display_network.m – Visualizes images or filters for autoencoders as a grid
  • pca_gen.m – Framework for whitening exercise
  • sampleIMAGESRAW.m – Returns 8×8 raw unwhitened patches

Softmax Regression |softmax_exercise.zip

  • checkNumericalGradient.m – Makes sure that computeNumericalGradient is implmented correctly
  • display_network.m – Visualizes images or filters for autoencoders as a grid
  • loadMNISTImages.m – Returns a matrix containing raw MNIST images
  • loadMNISTLabels.m – Returns a matrix containing MNIST labels
  • softmaxCost.m – Computes cost and gradient of cost function of softmax
  • softmaxTrain.m – Trains a softmax model with the given parameters
  • train.m – Framework for this exercise

from: http://ufldl.stanford.edu/wiki/index.php/MATLAB_Modules

人脸图像的几何归一化和灰度归一化

2017年10月13日

在对人脸表情进行识别时,人脸的归一化处理是至关重要的一环,它涉及到下一步处理的好坏。
人脸的归一化包括几何归一化和灰度归一化,几何归一化分两步:人脸校正和人脸裁剪。而灰度归
一化主要是增加图像的对比度,进行光照补偿。
1.几何归一化
几何归一化的目的主要是将表情子图像变换为统一的尺寸,有利于表情特征的提取。具体步骤如下:
(1)标定特征点,这里用[x,y] = ginput(3)函数来标定两眼和鼻子三个特征点。主要是用鼠标动手标
定,获取三个特征点的坐标值。
(2)根据左右两眼的坐标值旋转图像,以保证人脸方向的一致性。设两眼之间的距离为d,其中点为O。
(3)根据面部特征点和几何模型确定矩形特征区域,以O为基准,左右各剪切d,垂直方向各取0.5d和
1.5d的矩形区域进行裁剪。
(4)对表情子区域图像进行尺度变换为统一的尺寸,更有利于表情特征的提取。把截取的图像统一规格
为90*100的图像,实现图像的几何归一化。
面部几何模型如下图:
20130709122421625
2.灰度归一化
灰度归一化 主要是增加图像的亮度,使图像的细节更加清楚,以减弱光线和光照强度的影响。
这里用的是image=255*imadjust(C/255,[0.3;1],[0;1]); 用此函数进行光照补偿。
具体代码如下:

C= imread('Image001.jpg');
figure(1),imshow(C);
C=double(C);
image=255*imadjust(C/255,[0.3;1],[0;1]);
figure(2),imshow(image/255);
title('Lighting compensation');%光照补偿

[x,y] = ginput(3);    %%1 left eye, 2 right eye, 3 top of nose
cos = (x(2)-x(1))/sqrt((x(2)-x(1))^2+(y(2)-y(1))^2);
sin = (y(2)-y(1))/sqrt((x(2)-x(1))^2+(y(2)-y(1))^2);
mid_x = round((x(1)+x(2))/2);
mid_y = round((y(2)+y(1))/2);
d = round(sqrt((x(2)-x(1))^2+(y(2)-y(1))^2));
rotation = atan(sin./cos)*180/pi;
img = imrotate(image,rotation,'bilinear','crop'); 
figure(3), imshow(img);%人脸校正

[h,w] = size(img);
leftpad = mid_x-d;
if leftpad<1
   leftpad = 1;
end
toppad =mid_y - round(0.5*d);
if toppad<1
   toppad = 1;
 end
 rightpad = mid_x + d;
 if rightpad>w
    rightpad = w;
 end
 bottompad = mid_y + round(1.5*d);
 if bottompad>h
    bottompad = h;
 end   
 I1 =[];
 I2 =[];
 I1(:,:) = img(toppad:bottompad,leftpad:rightpad);
 I2(:,:) = imresize(I1,[90 100]); 
 figure(4),imshow(I2,[]);%人脸裁剪

机器学习模型需要对数据进行归一化

2017年10月11日

1)归一化后加快了梯度下降求最优解的速度;2)归一化有可能提高精度

1 归一化为什么能提高梯度下降法求解最优解的速度?
如下图所示,蓝色的圈圈图代表的是两个特征的等高线。其中左图两个特征X1和X2的区间相差非常大,X1区间是[0,2000],X2区间是[1,5],其所形成的等高线非常尖。当使用梯度下降法寻求最优解时,很有可能走“之字型”路线(垂直等高线走),从而导致需要迭代很多次才能收敛;而右图对两个原始特征进行了归一化,其对应的等高线显得很圆,在梯度下降进行求解时能较快的收敛。因此如果机器学习模型使用梯度下降法求最优解时,归一化往往非常有必要,否则很难收敛甚至不能收敛。
192105553858119
2 归一化有可能提高精度
一些分类器需要计算样本之间的距离(如欧氏距离),例如KNN。如果一个特征值域范围非常大,那么距离计算就主要取决于这个特征,从而与实际情况相悖(比如这时实际情况是值域范围小的特征更重要)。
3 归一化的类型
1)线性归一化
76512b142c1b7e27e8a7e7eb1fc11225
这种归一化方法比较适用在数值比较集中的情况。这种方法有个缺陷,如果max和min不稳定,很容易使得归一化结果不稳定,使得后续使用效果也不稳定。实际使用中可以用经验常量值来替代max和min。

2)标准差标准化
  经过处理的数据符合标准正态分布,即均值为0,标准差为1,其转化函数为:
z-score-2
  其中μ为所有样本数据的均值,σ为所有样本数据的标准差。

3)非线性归一化
经常用在数据分化比较大的场景,有些数值很大,有些很小。通过一些数学函数,将原始值进行映射。该方法包括 log、指数,正切等。需要根据数据分布的情况,决定非线性函数的曲线,比如log(V, 2)还是log(V, 10)等。
源码示例(这里包含了三种方法:前两种是编程实现,最后一种直接调用MATLAB函数来实现):

oriImage = imread('XXXX.jpg');
grayImage = rgb2gray(oriImage);
figure;
imshow(grayImage);

originalMinValue = double(min(min(grayImage)));
originalMaxValue = double(max(max(grayImage)));
originalRange = originalMaxValue - originalMinValue;

% Get a double image in the range 0 to +255
desiredMin = 0;
desiredMax = 255;
desiredRange = desiredMax - desiredMin;
dblImageS255 = desiredRange * (double(grayImage) - originalMinValue) / originalRange + desiredMin;

figure;
imshow(uint8(dblImageS255));

% Get a double image in the range 0 to +1
desiredMin = 0;
desiredMax = 1;
desiredRange = desiredMax - desiredMin;
dblImageS1 = desiredRange * (double(grayImage) - originalMinValue) / originalRange + desiredMin;

figure;
imshow(dblImageS1);

% Another way to normalazation, which only calls MATLAB function
img3 = mat2gray(oriImage);
figure;
imshow(img3);

MNIST数据库介绍及转换

2017年10月9日

数据库的一个子集。

MNIST数据库官方网址为:http://yann.lecun.com/exdb/mnist/ ,直接下载,train-images-idx3-ubyte.gz、train-labels-idx1-ubyte.gz等。下载四个文件,解压缩。解压缩后发现这些文件并不是标准的图像格式。这些图像数据都保存在二进制文件中。每个样本图像的宽高为28*28。

以下为将其转换成普通的jpg图像格式的代码:
Matlab

% Matlab_Read_t10k-images_idx3.m
% 用于读取MNIST数据集中t10k-images.idx3-ubyte文件并将其转换成bmp格式图片输出。
% 用法:运行程序,会弹出选择测试图片数据文件t10k-labels.idx1-ubyte路径的对话框和
% 选择保存测试图片路径的对话框,选择路径后程序自动运行完毕,期间进度条会显示处理进度。
% 图片以TestImage_00001.bmp~TestImage_10000.bmp的格式保存在指定路径,10000个文件占用空间39M。。
% 整个程序运行过程需几分钟时间。
% Written By DXY@HUST IPRAI
% 2009-2-22
clear all;
clc;
%读取训练图片数据文件
[FileName,PathName] = uigetfile('*.*','选择测试图片数据文件t10k-images.idx3-ubyte');
TrainFile = fullfile(PathName,FileName);
fid = fopen(TrainFile,'r'); %fopen()是最核心的函数,导入文件,‘r’代表读入
a = fread(fid,16,'uint8'); %这里需要说明的是,包的前十六位是说明信息,从上面提到的那个网页可以看到具体那一位代表什么意义。所以a变量提取出这些信息,并记录下来,方便后面的建立矩阵等动作。
MagicNum = ((a(1)*256+a(2))*256+a(3))*256+a(4);
ImageNum = ((a(5)*256+a(6))*256+a(7))*256+a(8);
ImageRow = ((a(9)*256+a(10))*256+a(11))*256+a(12);
ImageCol = ((a(13)*256+a(14))*256+a(15))*256+a(16);
%从上面提到的网页可以理解这四句
if ((MagicNum~=2051)||(ImageNum~=10000))
    error('不是 MNIST t10k-images.idx3-ubyte 文件!');
    fclose(fid);    
    return;    
end %排除选择错误的文件。
savedirectory = uigetdir('','选择测试图片路径:');
h_w = waitbar(0,'请稍候,处理中>>');
for i=1:ImageNum
    b = fread(fid,ImageRow*ImageCol,'uint8');   %fread()也是核心的函数之一,b记录下了一副图的数据串。注意这里还是个串,是看不出任何端倪的。
    c = reshape(b,[ImageRow ImageCol]); %亮点来了,reshape重新构成矩阵,终于把串转化过来了。众所周知图片就是矩阵,这里reshape出来的灰度矩阵就是该手写数字的矩阵了。
    d = c'; %转置一下,因为c的数字是横着的。。。
    e = 255-d; %根据灰度理论,0是黑色,255是白色,为了弄成白底黑字就加入了e
    e = uint8(e);
    savepath = fullfile(savedirectory,['TestImage_' num2str(i,'d') '.bmp']);
    imwrite(e,savepath,'bmp'); %最后用imwrite写出图片
    waitbar(i/ImageNum);
end
fclose(fid);
close(h_w);

CPP

#include "funset.hpp"
#include <iostream>
#include <fstream>
#include <vector>
#include <opencv2/opencv.hpp>

static int ReverseInt(int i)
{
	unsigned char ch1, ch2, ch3, ch4;
	ch1 = i & 255;
	ch2 = (i >> 8) & 255;
	ch3 = (i >> 16) & 255;
	ch4 = (i >> 24) & 255;
	return((int)ch1 << 24) + ((int)ch2 << 16) + ((int)ch3 << 8) + ch4;
}

static void read_Mnist(std::string filename, std::vector<cv::Mat> &vec)
{
	std::ifstream file(filename, std::ios::binary);
	if (file.is_open()) {
		int magic_number = 0;
		int number_of_images = 0;
		int n_rows = 0;
		int n_cols = 0;
		file.read((char*)&magic_number, sizeof(magic_number));
		magic_number = ReverseInt(magic_number);
		file.read((char*)&number_of_images, sizeof(number_of_images));
		number_of_images = ReverseInt(number_of_images);
		file.read((char*)&n_rows, sizeof(n_rows));
		n_rows = ReverseInt(n_rows);
		file.read((char*)&n_cols, sizeof(n_cols));
		n_cols = ReverseInt(n_cols);

		for (int i = 0; i < number_of_images; ++i) {
			cv::Mat tp = cv::Mat::zeros(n_rows, n_cols, CV_8UC1);
			for (int r = 0; r < n_rows; ++r) {
				for (int c = 0; c < n_cols; ++c) {
					unsigned char temp = 0;
					file.read((char*)&temp, sizeof(temp));
					tp.at<uchar>(r, c) = (int)temp;
				}
			}
			vec.push_back(tp);
		}
	}
}

static void read_Mnist_Label(std::string filename, std::vector<int> &vec)
{
	std::ifstream file(filename, std::ios::binary);
	if (file.is_open()) {
		int magic_number = 0;
		int number_of_images = 0;
		int n_rows = 0;
		int n_cols = 0;
		file.read((char*)&magic_number, sizeof(magic_number));
		magic_number = ReverseInt(magic_number);
		file.read((char*)&number_of_images, sizeof(number_of_images));
		number_of_images = ReverseInt(number_of_images);

		for (int i = 0; i < number_of_images; ++i) {
			unsigned char temp = 0;
			file.read((char*)&temp, sizeof(temp));
			vec[i] = (int)temp;
		}
	}
}

static std::string GetImageName(int number, int arr[])
{
	std::string str1, str2;

	for (int i = 0; i < 10; i++) {
		if (number == i) {
			arr[i]++;
			str1 = std::to_string(arr[i]);

			if (arr[i] < 10) {
				str1 = "0000" + str1;
			} else if (arr[i] < 100) {
				str1 = "000" + str1;
			} else if (arr[i] < 1000) {
				str1 = "00" + str1;
			} else if (arr[i] < 10000) {
				str1 = "0" + str1;
			}

			break;
		}
	}

	str2 = std::to_string(number) + "_" + str1;

	return str2;
}

int MNISTtoImage()
{
	// reference: http://eric-yuan.me/cpp-read-mnist/
	// test images and test labels
	// read MNIST image into OpenCV Mat vector
	std::string filename_test_images = "E:/GitCode/NN_Test/data/database/MNIST/t10k-images.idx3-ubyte";
	int number_of_test_images = 10000;
	std::vector<cv::Mat> vec_test_images;

	read_Mnist(filename_test_images, vec_test_images);

	// read MNIST label into int vector
	std::string filename_test_labels = "E:/GitCode/NN_Test/data/database/MNIST/t10k-labels.idx1-ubyte";
	std::vector<int> vec_test_labels(number_of_test_images);

	read_Mnist_Label(filename_test_labels, vec_test_labels);

	if (vec_test_images.size() != vec_test_labels.size()) {
		std::cout << "parse MNIST test file error" << std::endl;
		return -1;
	}

	// save test images
	int count_digits[10];
	std::fill(&count_digits[0], &count_digits[0] + 10, 0);

	std::string save_test_images_path = "E:/GitCode/NN_Test/data/tmp/MNIST/test_images/";

	for (int i = 0; i < vec_test_images.size(); i++) {
		int number = vec_test_labels[i];
		std::string image_name = GetImageName(number, count_digits);
		image_name = save_test_images_path + image_name + ".jpg";

		cv::imwrite(image_name, vec_test_images[i]);
	}

	// train images and train labels
	// read MNIST image into OpenCV Mat vector
	std::string filename_train_images = "E:/GitCode/NN_Test/data/database/MNIST/train-images.idx3-ubyte";
	int number_of_train_images = 60000;
	std::vector<cv::Mat> vec_train_images;

	read_Mnist(filename_train_images, vec_train_images);

	// read MNIST label into int vector
	std::string filename_train_labels = "E:/GitCode/NN_Test/data/database/MNIST/train-labels.idx1-ubyte";
	std::vector<int> vec_train_labels(number_of_train_images);

	read_Mnist_Label(filename_train_labels, vec_train_labels);

	if (vec_train_images.size() != vec_train_labels.size()) {
		std::cout << "parse MNIST train file error" << std::endl;
		return -1;
	}

	// save train images
	std::fill(&count_digits[0], &count_digits[0] + 10, 0);

	std::string save_train_images_path = "E:/GitCode/NN_Test/data/tmp/MNIST/train_images/";

	for (int i = 0; i < vec_train_images.size(); i++) {
		int number = vec_train_labels[i];
		std::string image_name = GetImageName(number, count_digits);
		image_name = save_train_images_path + image_name + ".jpg";

		cv::imwrite(image_name, vec_train_images[i]);
	}

	// save big imags
	std::string images_path = "E:/GitCode/NN_Test/data/tmp/MNIST/train_images/";
	int width = 28 * 20;
	int height = 28 * 10;
	cv::Mat dst(height, width, CV_8UC1);

	for (int i = 0; i < 10; i++) {
		for (int j = 1; j <= 20; j++) {
			int x = (j-1) * 28;
			int y = i * 28;
			cv::Mat part = dst(cv::Rect(x, y, 28, 28));

			std::string str = std::to_string(j);
			if (j < 10)
				str = "0000" + str;
			else
				str = "000" + str;

			str = std::to_string(i) + "_" + str + ".jpg";
			std::string input_image = images_path + str;

			cv::Mat src = cv::imread(input_image, 0);
			if (src.empty()) {
				fprintf(stderr, "read image error: %s\n", input_image.c_str());
				return -1;
			}

			src.copyTo(part);
		}
	}

	std::string output_image = images_path + "result.png";
	cv::imwrite(output_image, dst);

	return 0;
}

LDA算法入门

2017年9月29日
参考
Introduction to LDA
Linear Discriminant Analysis – A Brief Tutorial
http://www.aiaccess.net/English/Glossaries/GlosMod/e_gm_fisher_discriminant.htm
线性判别分析(Linear Discriminant Analysis, LDA)算法分析

一. LDA算法概述:

线性判别式分析(Linear Discriminant Analysis, LDA),也叫做Fisher线性判别(Fisher Linear Discriminant ,FLD),是模式识别的经典算法,它是在1996年由Belhumeur引入模式识别和人工智能领域的。性鉴别分析的基本思想是将高维的模式样本投影到最佳鉴别矢量空间,以达到抽取分类信息和压缩特征空间维数的效果,投影后保证模式样本在新的子空间有最大的类间距离和最小的类内距离,即模式在该空间中有最佳的可分离性。因此,它是一种有效的特征抽取方法。使用这种方法能够使投影后模式样本的类间散布矩阵最大,并且同时类内散布矩阵最小。就是说,它能够保证投影后模式样本在新的空间中有最小的类内距离和最大的类间距离,即模式在该空间中有最佳的可分离性。

下载详细介绍文件文档  以及 文档和代码

» 阅读更多: LDA算法入门