3.2 一般矩阵运算
前面介绍了创建矩阵的基本方法,现在我们来看一些常用的矩阵运算,包括线性代数运算、矩阵索引和矩阵元素筛选。
3.2.1 线性代数运算
你可以对矩阵进行各种线性代数运算,比如矩阵相乘、矩阵数量乘法和矩阵加法。针对以前例子中的y,以下为这三种运算的实例:
关于矩阵线性代数运算的更多细节请参见8.4节。
3.2.2 矩阵索引
2.4.2节中的向量运算同样适用于矩阵,例如:
这里我们提取了矩阵z中第2、3列的所有元素组成了一个子矩阵。
下面的例子提取的是矩阵的行:
还可以对一个矩阵的子矩阵进行赋值:
这里给y的第1、3行赋了新的值。
下面是另一个给子矩阵赋值的例子:
向量的负值索引用来排除掉某些元素,这种操作同样适用于矩阵。
命令第二行的意思是,取出矩阵y中除第二行外的所有行。
3.2.3 扩展案例:图像操作
图像文件本质上就是矩阵,因为像素点也都是按行和列排列的。一张灰度图像会把每一个像素点的亮度存储为矩阵的一个元素。例如,图像中位于第28行、第88列的像素点的灰度值就存储在矩阵的第28行、第88列。如果图像是彩色的话,就需要三个矩阵来存储,分别记录红、黄、蓝的强度值。这里我们只以灰度图像为例。
以一张拉什莫尔山国家纪念公园的照片为例,我们用pixmap包读取图像。(附录B介绍了下载和安装R包的方法。)
读取文件mtrush1.pgm,并返回一个pixmap类的对象。再把这个对象在R里画出来,如图3-1所示。
我们来看一下pixmap类的组成内容:
pixmap类属于S4类型,它用符号@来访问各个组件,而不是用符号$。S3和S4的问题将在第9章讨论,不过这里最关键的是灰度矩阵mtrush1@grey。在这个例子里,这个矩阵有194行和259列。
这个例子的灰度是用0.0到1.0之间的数值表示的,0.0表示黑色,1.0表示白色,中间值代表不同程度的灰色。例如,图片第28行、第88列的点非常亮。
为了演示矩阵操作,我们将罗斯福总统的脸覆盖上(对不起了,总统先生,我对你没有个人恩怨)。使用R中locator()函数来找到需要盖住部分的行号和列号。当调用这个函数时,R会等待用户点击图片上的任意点,然后返回该点的准确坐标。用这种方式可以找到罗斯福总统的脸对应区域是从84行到163行,从135列到177列。(注意,pixmap对象的行编号和locator()是相反的:pixmap是自上而下递增的,而locator()恰恰相反。)于是,为了把这部分图像覆盖上,我们把这一区域的所有像素点取值设为1.0。
效果如图3-2所示。
如果只是想将罗斯福总统的脸模糊掉,而不是去掉,可以在图片上增加随机噪声。代码如下:
如注释所显示的,我们生成随机噪声,然后把目标像素点矩阵和噪声进行加权平均。参数q控制噪声的权重,q值越大就越模糊。随机噪声取自U(0,1),即区间(0,1)上的均匀分布。注意,下面是矩阵操作:
效果如图3-3所示。
3.2.4 矩阵元素筛选
和向量一样,矩阵也可以做筛选。但是需要注意一下语法上的不同。首先来看一个简单的例子:
下面来分析一下这条命令,正如我们在第2章所讲那样:
为了提高性能,需要注意的是,计算j时使用的是完全向量化运算。这是因为:
x[,2]是向量。
运算符>=用于比较两个向量。
数值3被自动重复,变成一个由数值3组成的向量。
还需要注意,虽然这个例子中j是通过x定义且用于提取x中的元素,但事实上筛选规则可以基于除被筛选变量之外的变量。以下例子同样使用上述的x:
表达式 z 2 ==1用于判断z的每个元素是否是奇数,返回的结果是(TRUE, FALSE, TRUE)。因此我们提取出了x的第一、三行。
再看另外一个例子:
这里也是同样的做法,但是使用了稍微复杂的条件来提取矩阵的行。(用类似的方式可以提取列或子矩阵。)首先表达式m[,1] > 1把m的第一列的每个元素同1进行比较,返回(FALSE,TRUE,TRUE)。类似的是,第二个表达式m[,2] > 5返回(FALSE,FALSE,TRUE)。对(FALSE,TRUE,TRUE)和(FALSE,TRUE,TRUE)进行逻辑“与”运算,得到(FALSE,FALSE,TRUE)。以运算后的向量为行索引,得到矩阵m的第三行。
注意,这里需要使用 & 运算符,而不是 &&,前者是向量的逻辑“与”运算,后者是用于if语句的标量逻辑“与”运算。想了解更多可以参见7.2节中的运算符列表。
细心的读者可能注意到前面的例子中有一些反常的地方。我们要得到的是一个1行2列的子矩阵,可是返回的却是一个由两个元素组成的向量。元素是正确的,可是数据类型却不对。如果把返回值输入到其他的矩阵函数里可能会导致错误。解决办法是通过设定参数drop,告诉R让返回结果保持数据的二维属性。3.6节会介绍如何避免意外降维,并详细讲解drop的用法。
由于矩阵也是向量,所以向量的运算也适用于矩阵。例如:
这行命令说明,从向量索引的角度来看,m的第1、3、5、6个元素大于2。第5个元素指的是m第2行第2列的元素,即10,的确比2大。
3.2.5 扩展案例:生成协方差矩阵
这个例子演示row()和col()函数,这两个函数的参数都是矩阵。例如对于矩阵a,row(a[2,8])返回a中对应元素的行号,即2。然而,显然我们知道row(a[2,8])在第2行,那么这个函数有什么用途呢?
先来看一个例子。如果我们想要用MASS库的mvrnorm()函数生成多元正态分布的随机数,需要指定协方差矩阵。协方差矩阵是对称的,就是说矩阵中第1行第2列的元素等于第2行第1列的元素。
假如我们在处理一个n元正态分布,它的协方差矩阵有n行n列。我们希望这n个随机变量方差都为1,每两个变量间的相关性都是rho。例如当n=3,rho=0.2时,需要的矩阵就是:
下面就是生成这种矩阵的代码:
下面来解释它的运算过程。首先col()返回的是矩阵元素的列号,和row()功能相似。然后第三行的row(m)返回一个由整数组成的矩阵,每个整数代表矩阵m中对应元素的行号。例如:
因此表达式row(m) == col(m)返回一个由FALSE和TRUE组成的矩阵,矩阵的对角线位置是TRUE,其他位置是FALSE。再次提醒一下读者,本例中的二目运算符"=="是一个函数,row()和col()也是函数。因此,下面这个表达式:
会作用到矩阵m的每一个元素上,最后返回一个由FALSE和TRUE组成的矩阵,其行数和列数与m的相同。表达式ifelse()同样是一个函数调用。
在本例中,刚才提到的由TRUE和FALSE组成的矩阵传递给ifelse()函数做参数,返回一个矩阵,对角线元素为1,其他位置为0。