《R语言编程艺术》——3.2　一般矩阵运算

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时，需要的矩阵就是：