呆在家无聊？何不抓住这个机会好好学习！

本公众号一向坚持的理念是数据分析工具要从基础开始学习，按部就班，才能深入理解并准确利用这些工具。鼠年第一篇原创推送比较长，将从基础的线性代数开始。线性代数大家都学过，但可能因为联系不到实用情况，都还给了曾经的老师。线性代数是数理统计尤其是各种排序分析的基础，今天我将以全新的角度基于R语言介绍线性代数，并手动完成PCA分析，从而强化关于线性代数和实际数据分析的联系。

就像把标量x变成2x一样，我们经常需要把向量(x, y)变成诸如(2x+y, x–3y)的新向量，也即线性变换。那么可以定义矩阵乘法，用于表示一切线性变换，对一个行向量进行线性变换就是右乘变换矩阵，对一个列向量进行线性变换就是左乘变换矩阵。从几何上看，把平面上的每个点(x, y)都变到(2x+y, x–3y)的位置上去，效果就相当于对这个平面进行了一个“线性的拉扯”。

矩阵的乘法，其实就是多个线性变换叠加的效果，它显然满足结合律，但不满足交换律（满足交换律的东西要求运算符两边具有等质性，而线性变换左边为变换对象，右边为变换方式，显然不满足等质性）。主对角线全是1其余元素为0的矩阵所对应的线性变换其实就是不变的意思，因此它叫做单位矩阵。矩阵A乘以矩阵B得单位矩阵，就是做完线性变换A后再做一次线性变换B就又变回原对象，因此矩阵B是矩阵A的逆矩阵。一个矩阵是可逆的就意味着经过这个矩阵变换的图形可以还原。

矩阵的秩为线性变换的维度，方阵对应的行列式的绝对值是每个单位正方形在经过该方阵变换之后的面积，或者任意图形经过该方阵变换之后面积变化的倍数（伸缩因子），行列式值为负改变基向量的相对位置。因此，单位矩阵的行列式值为1，某行全为0或秩小于n的n阶矩阵的行列式值为0（因为某些维度会被删除，线性变换会导致降维），而且|A·B|显然等于|A|·|B|。行列式值为0，那么对应的矩阵当然不可逆，因为这样的线性变换会导致降维，无法进行还原。因此可逆的矩阵一定是方阵，且其对应的行列式值不为0。

行列式只表征了矩阵变换前后面积的变化，对于方阵来说，得知其变换中运动的“方向”与“距离”十分重要。矩阵变换的“方向”为矩阵的特征向量，其“距离”为矩阵的特征值，特征值与特征向量描绘了矩阵变换的矢量特征。因此，线性代数即是研究向量、向量空间（线性空间）、线性变换和有限维线性方程组的数学分支。由于科学研究所获得的数据往往可以用不同向量构成的数据集来表示，线性对数对于数据降维、去冗余等有十分重要的应用。

任何一个数据向量是不可能孤立存在的，必须基于一定的坐标系，只不过通常默认的是单位矩阵所代表的规范正交坐标系。

目前人类认识最全面的系统就是线性系统，线性代数是区别于微积分这个系统的另外一个描述世界的最有力武器。线性代数在系统力学、统计学、运筹学、机器学习、数值计算、图像处理、工程设计、信号处理、自动化控制等领域均有十分重要的应用，几乎渗透到计算机应用的各个方面。线性代数是现代数学的语言，其不仅仅是一个工具，更是一种思维，著名数学软件Matlab的全称就是Matrix laboratory（矩阵实验室）。

矩阵与行列式

向量、矩阵与行列式是线性代数研究的基本对象，注意这里的矩阵为数学概念，与R语言中的矩阵不能等同，但是数学中的矩阵可以利用R中的矩阵来存储，例如在R中可以用函数matrix()来创建一个矩阵：

当然，也可以使用其他任何来源的数据创建与储存矩阵，这里主要讨论数学中矩阵与行列式的运算及其在R中的实现。

⑴向量的运算

向量运算是矩阵运算的基础，向量的加减法为对应元素的加减。一个行向量点乘以一个同维列向量称作向量的内积，又叫作点积，结果是一个数；一个列向量叉乘以一个行向量称作向量的外积，结果是一个矩阵。假如a=(a1,a2)和b=(b1,b2)为两个列向量，那么点乘与叉乘的区别如下所示：

点乘可以理解为降维运算，在R中的符号位%*%，也可以使用crossprod()函数；叉乘为升维运算，在R中可以使用tcrossprod()函数，如下所示：

内积为一个向量在另一个向量方向上的投影，因此在解析几何中有如下等式：

a·b=|a||b|cosθ

其中θ为两个向量的夹角，不难想象假如a·b=0，那么θ为90˚，a与b正交。

⑵矩阵的运算

具有m行n列的矩阵称为m×n矩阵，共具有m×n个元素；行和列数均为n的称为n阶矩阵或n阶方阵。只有一行的矩阵为行向量，只有一列的矩阵为列向量，行数和列数均相等的矩阵称为同型矩阵。不在对角线上的元素均为0的方阵，则称为对角矩阵，对角线上元素均为1的对角矩阵为单位矩阵，记作E。

①矩阵与数的四则运算

矩阵与数的四则运算为每个元素与数的四则运算，其结果为原来矩阵的同型矩阵，为矩阵的线性运算，如下所示：

②矩阵与矩阵的加/减法

矩阵与矩阵的加/减法为对应元素的两两相加/减（必须是同型矩阵之间），满足交换律和结合律，同为矩阵的线性变换，如下所示：

③矩阵与矩阵的乘法

设A为m×p的矩阵，B为p×n的矩阵，那么称m×n的矩阵C为矩阵A与B的乘积，记作C=AB，其中矩阵C中的第i行第j列元素可以表示为：

实例如下所示：

矩阵与矩阵相乘不满足结合律，但是满足交换律和分配律，在R中可使用%*%符号来计算，如下所示：

矩阵相乘的Hadamard乘积定义为矩阵每个对应元素的乘积（必须是两个同型矩阵之间），在R中使用*符号来计算：

④对角有关的运算

把矩阵的行换成列称为矩阵的转置，如果矩阵A的转置矩阵等于本身也即AT=A，那么称之为对称矩阵，对角矩阵一定为对称阵。在R中矩阵转置可以使用t()函数，diag(v)表示以向量v的元素为对角线元素的对角阵，当M是一个矩阵时，则diag(M)表示的是取M对角线上的元素构造向量，如下所示：

在R中，我们可以很方便的取到一个矩阵的上、下三角部分的元素，函数lower.tri()和函数upper.tri()提供了有效的方法。例如lower.tri(x, diag =FALSE)，那么函数将返回一个逻辑值矩阵，其中下三角部分为真，上三角部分为假，选项diag为真时包含对角元素，为假时不包含对角元素。upper.tri()则与之相反，取矩阵上三角部分，具体如下所示：

⑤与维数有关

在R中很容易得到一个矩阵的维数（指矩阵的行数和列数），函数dim()将返回一个矩阵的维数，此外nrow()和ncol()分别返回行数和列数，row()和col()则返回矩阵每个元素的行数与列数坐标，如下所示：

⑶行列式的运算

由n阶方阵A的元素构成的行列式，称为方阵A的行列式，记作|A|或者detA，在R中函数det()是求矩阵行列式的值，如下所示：

假如存在矩阵A*使得AA*=A*A=|A|E，那么A*为A的伴随矩阵，而如果存在矩阵B使得AB=BA=E，那么称B为A的逆矩阵（即A-1），逆矩阵和伴随矩阵的关系如下所示：

可以看出若矩阵可逆那么其行列式值一定不能为0。在R中矩阵求逆可用函数solve()，应用solve(a, b)运算结果是解线性方程组ax=b，若b缺省，则系统默认为单位矩阵，因此可用其进行矩阵求逆，例如：

线性变换

线性变换可以用矩阵表示，那么如何描述线性变换的特征，需要用到矩阵的一些属性。我们已经知道矩阵的行列式为线性变换的伸缩因子，其他矩阵属性还有矩阵的秩、特征值、特征向量等。

⑴矩阵的秩

矩阵可以理解为一组线性方程的系数矩阵，线性方程组为线性变换的具体形式，任何一组如下所示的有限数目线性方程：

都可以表示成Ax=b的形式，其中A为系数矩阵，x为未知数向量或解向量，b为常数向量。根据矩阵转置运算规则（Ax）T=xTAT=bT，因此方程组的实际含义为求解能经过A变换得到b的列向量或者AT变换得到bT的行向量。克拉默法则告诉我们，如果n=m（也即A为方阵）且|A|不为0的话那么线性方程组一定有解。然而现实中n=m的例子并不多见，但是我们知道，当方程数目过多的时候不同的方程之间一点存在线性相关性，可以通过方程式间的初等运算（比如数乘加减等）消除掉冗余的方程，从而获得最精简的方程组来求解，表现在矩阵上就是矩阵B=（A，b）的初等变换。在一个m×n矩阵中任取k行k列元素构成的行列式称为矩阵的k阶子式。经过初等运算后最精简的线性方程的个数也即经过初等变换后矩阵B的最高阶非零子式的阶数称之为矩阵B的秩（rank），记为R(B)，初等变换也即保持秩不变的线性变换。根据A和B的秩的大小可以判断是否存在列向量可以经过A变换得到b。一个变换矩阵的秩可以理解为图像经过该矩阵变换之后的维度。因此如果B的秩大于A，也即结果向量b的维度高于变换矩阵A，方程组一定无解。

对于n阶矩阵A，如果R(A)=n，则为满秩矩阵，且有|A|不为0，矩阵可逆；如果R(A)<n，则为降秩矩阵，且有|A|=0，矩阵不可逆（称为奇异矩阵）。矩阵不可逆则其线性变换不可逆。

⑵向量空间

空间的一个基本属性就是容纳符合规则的运动（变换），线性代数所研究的空间为向量空间，在向量空间可以进行线性变换，矩阵的本质就是描述运动，不过这里的运动与微积分描述的连续运动不同，是瞬间的没有过程的跃迁。向量中元素的个数为维数，多个同维向量可以组成向量组。假如给定向量组A:a1,a2, … am，如果存在不全为0的数k1,k2, … km，使得：

k1a1+k2a2+…+kmam=0

则称向量组A线性相关，否则为线性无关。向量组A内的最大无关向量组称之为向量组的秩，向量组内的向量均可用最大无关向量组内的向量进行线性表示。向量组A的秩等于矩阵A的秩，那么就有R(A)≤n，假如R(A)<m，A肯定线性相关。在线性空间中，全体n维向量所组成的集合，称为n维向量空间，记为Rn，假如含有n个线性无关的向量的向量组A:a1, a2, … am，使得Rn内的所有向量均可使用这n个向量线性表示，则称A为Rn的一个基，假如a1,a2, … am两两正交，那么则称A为Rn的正交基，假如a1, a2, … am都是单位向量，那么则称A为Rn的规范正交基，把一个普通基转换为规范正交基的过程称为规范正交化。

在向量空间，变换的对象是向量，向量可以用这个向量空间的一个基（也即坐标系）来描述，一个向量空间有无数个向量，也有无数个基，每个基实际上都是一个非奇异矩阵，由此看来对一个向量的矩阵变换（线性变换）实际上就是坐标系的变换。这时候我们回到最开头的例子，向量(x,y)实际上默认的是

这个规范正交基描述的，对其使用

进行变换成(2x+y,x–3y)，可以理解为固定坐标系对向量进行了变换，也可以理解为固定向量对坐标系进行了变换，这与运动是相对的思想是符合的。

推广而来，对于一个列向量a进行矩阵M变换则为Ma，实际上是用M坐标系去描述a，而a实际上默认是使用的坐标系E进行描述的，换言之在向量空间没有单独存在的向量，我们看到的向量a实际上是Ea，E是最特殊的一个规范正交基，其他基均可通过初等变换转换成E，也即其他基均可用E进行描述。这是可以理解的，在空间中向量不是单独存在的事物，向量本质上就是在坐标系中不同坐标的投影值，矩阵变换其实就是坐标系的变换。因此，Ma=b实际上就是MEa=Eb，也即将E中的a转换到M中，就是E中的b，也即最终的结果都是要在E中描述出来。矩阵的乘法也即不同线性变换的叠加，或者使用一个坐标系去描述另一个坐标系。

假如n个n维向量是n维向量空间的一个规范正交基，那么向量组对应的n阶矩阵A称之为正交矩阵，正交阵满足下面关系：

ATA=E(即A-1=AT)

正交矩阵对应的为正交变换，且|A|=1或者-1，也即正交变换能保持图形的形状大小不变，正交变换实际上为坐标系的保形旋转。

⑶特征值与特征向量

设有n阶矩阵A，如果存在数λ和n维非零列向量x使得Ax=λx，那么称λ为矩阵A的特征值，x称为矩阵A对应与特征值λ的特征向量。Ax=λx实际上就是Ax=λEx，也就是说向量x使用A来描述时仅仅是长度的变化而无方向的旋转，唯一可能的解释就是矩阵A的变换的运动方向就是沿着x的方向。

但是我们不能急于下结论，再回到开头的例子（如下图所示），可以看到，矩阵变换可能会导致图形的形状与大小的改变，不同点的运动方向并不一样：

矩阵A的特征值可通过特征方程|A-λE|=0进行求解，例如上图矩阵可求得其特征根为

其约为2.19和-3.19，其对应的特征向量可分别表示为k（2.5+0.5，1）和k（2.5-0.5，1），也即约为k（5.19，1）和k（-0.19，1），其中k为不等0的常数。特征向量不唯一，但是同一特征值对应的特征向量其方向相同，位置关系如上图所示，可以看出两个特征向量是正交的。在向量的矩阵变换中，不同的向量变换的方向、距离不一样，但是矩阵特征值λ对应的特征向量其变换方向不变，仅进行比例为λ的长度伸缩。

特征向量对于矩阵变换方向的指示作用在于：任何一个向量在进行变换后其在特征值λ对应的特征向量方向上的投影均缩放为原来的λ倍。特征值λ为正，投影方向不发生变化；特征值λ为负，投影的方向则发生变化（也即从原点一侧到另一侧。对于实数矩阵A来说，不同实特征值的个数等于其矩阵的秩，且不同特征值对应的特征向量之间相互正交，满秩矩阵全部特征值的乘积等于其行列式|A|的值，满秩矩阵的特征向量矩阵为正交矩阵。

⑷相似矩阵

当在一个向量空间使用一个坐标系（也即一组线性无关的向量）来描述一个向量时，实际上使用的是这个向量在其他向量上的投影。我们可以使用任意坐标系来进行描述，然而很多坐标系其向量组可能不是两两正交的，显然两两正交的向量组更适合做坐标系。假如把矩阵变换看成坐标系变换，也即使用矩阵的坐标系来描述向量，那么求特征向量实际上是一个正交化过程，使用特征向量上的投影也即特征值来描述所对应的坐标系。根据特征值与特征向量的定义：Ap=λp，对于非奇异矩阵，假设特征向量构成的特征矩阵为P，那么上式有AP=PΛ，也即P-1AP=Λ，其中Λ为对角线上是特征值的对角矩阵。假如矩阵B与A具有相同的特征值，其对应的特征向量矩阵为Q，那么有：

Q-1BQ=Λ=P-1AP

B=QP-1APQ-1

又由于(PQ-1)-1=QP-1，因此上式可以一般化为B=M-1AM，也即尽管A和B具有不同的特征矩阵（正交化的坐标系），但是A和B在各自特征向量上的投影也即特征值相同，而这两个正交化的特征向量坐标系是可以通过简单旋转来转换的（因为P、Q均为正交矩阵，也即正交转换），我们称B为A的相似矩阵，M为A到B的相似变换矩阵，此外矩阵A与由其特征值组成的对角矩阵Λ也相似，上式P-1AP=Λ将A转换为对角矩阵的线性变换过程称为矩阵对角化。

现在我们返回A=PΛP-1这个式子，由于P为正交矩阵，进行的是坐标系旋转变换，Λ为对角阵，进行的是拉伸变换，因此上式是对A矩阵变换的分解，称为特征分解（Eigendecomposition）、特征值分解（Eigen ValueDecomposition，EVD）或谱分解（Spectraldecomposition）。因此，相似矩阵可以理解为拥有相同的拉伸构型但有不同的旋转。

主成分分析详解

主成分分析（PrincipalComponent Analysis，PCA），是一种数理统计方法。通过正交变换将一组可能存在相关性的变量转换为一组线性不相关的变量，实现数据降维，转换后的这组变量叫主成分。

以下面水质监测数据为例进行分析：

water_quality=read.csv("water_quality.csv", row.names=1)
#数据下载链接：https://pan.baidu.com/s/1reAwppb6XFyYY0ld7kDHoQ 提取码：lzic

此数据中共有10个样点5个水质指标的检测结果，难以将其直接进行排序展示，需要进行降维，降维的目的就是“降噪”和“去冗余”。在阐述PCA的原理之前，首先需要明白的是，一个维度（也即变量）在不同样本间的方差越大，则对区分不同样本的贡献越大，也即对数据的解释能力越大。

假设数据中某个主要的维度A，它能代表原始数据，是“我们真正想听到的东西”，它本身含有的“能量”(即该维度的方差)本来应该是很大的，但由于它与其他维度存在一定的相关性，受到这些个相关维度的干扰，使得它的能量被削弱了，我们就希望通过PCA处理后，使维度A与其他维度的相关性尽可能减弱，进而恢复维度A应有的解释能力。

冗余也就是多余的意思，就是有它没它都一样。同样，假如数据中有些个维度，在所有的样本上变化不明显(极端情况：在所有的样本中该维度都等于同一个数)，也就是说该维度上的方差接近于零，那么显然它对区分不同的样本丝毫起不到任何作用，这个维度即是冗余的。

要进行降维，首先需要获得数据的协方差矩阵，协方差矩阵度量的是维度与维度之间的关系，而非样本与样本之间。协方差矩阵的主对角线上的元素是各个维度上的方差，其他元素是两两维度间的协方差(衡量相关性)。那么“降噪”，就是让保留下的不同维度间的相关性尽可能小，也就是说让协方差矩阵中非对角线元素都基本为零，也即矩阵对角化过程。而对角化后得到的矩阵，其对角线上是协方差矩阵的特征值，同时也是原始维度经过线性变换之后新维度上的方差，因此“去冗余”，则需要去掉特征值小也即方差小或者说对样本解释量小的维度，筛选高特征值对应的维度也即主成分。因此，PCA的原理即求解和筛选协方差矩阵特征值与特征向量。

假设原始数据矩阵为A（本例中列为所研究的维度，行为样本），将每个维度标准化后矩阵为C，协方差矩阵为S，那么S=CTC，在R中方法如下所示：

zdata=scale(water_quality, center=T, scale=T)
cov=cov(zdata)

可以看到对角线上都是1，这主要是因为对不同变量均进行了z-score标准化（均值为0，方差为1），这是一种妥协策略，因为不同变量的量纲不同，很难判断不同变量的初始解释能力，因此赋予每一个变量相等的方差。根据数据类型不同也可尝试其他标准化方法，例如没有离群点时可尝试Min-Max标准化等。接下来求协方差矩阵的特征根与特征向量：

cov.eigen=eigen(cov)

其中values为特征值，vectors为特征向量构成的特征矩阵设为P（正交矩阵），那么我们可以使用P-1SP对角化来检验这个过程：

t(cov.eigen$vectors)%*%cov%*%cov.eigen$vectors

可以看到主对角线即为特征值，其余位置元素接近0。对角化后的矩阵为变换后新维度的协方差矩阵，再根据下式：

Λ=P-1SP=PTSP=PTCTCP=(CP)TCP

也即Λ为CP的协方差矩阵，现进行计算如下：

newdim=zdata%*%cov.eigen$vectors

以上为样本在新维度上的坐标，由于后面三个特征值远小于前两个，可以选取前两个维度作为主成分绘制排序图，用图上样品点之间的欧氏距离表征样品之间的差异，如下所示：

plot(newdim[,1:2],xlab="PC1",ylab="PC2",pch=1:10)

其中CP为对C的行向量进行线性变换，那么新维度是列也即水质参数的线性组合，其系数也即P的列向量，我们可以对比哪一个水质参数对主成分的影响大以及如何影响，如下所示：

coef=data.frame(PC1=cov.eigen$vectors[,1],PC2=cov.eigen$vectors[,2],row.names=colnames(zdata))

主成份1与叶绿素a、总氮、总磷、COD、浮游植物量均正相关，系数为正且相差不大，而这几个指标越大均表示富营养化越严重，因此主成分1表示的是富营养化程度。主成分2前四个指标系数为负，浮游植物系数为正，且最后一个指标系数远大于前四个，可以看出主成分2主要表示浮游植物的生物量的增加以及浮游植物对营养元素尤其是氮、磷的消耗，也即在水体富营养化后期。

PCA是基于维度（也即变量）之间的协方差矩阵进行分析，实际上PCA只是进行了维度的正交化并给出正交化后每个维度的贡献（特征值），正交化的维度也即主成分其个数等于原来数据矩阵的秩，之后根据新维度方差贡献的大小而忽略贡献率小的坐标、保留方差贡献率大的坐标，从而实现了数据的降维。

此PCA实例来自我研究生数学老师袁卓建老师，在此像袁老师致敬！基于特征值分解的PCA分析是很多降维排序分析的基础，例如主坐标分析（Principal Coordinate Analysis，PCoA）、冗余分析(Redundancy analysis，RDA)等。