@evilking
2017-10-15T10:38:24.000000Z
字数 11667
阅读 2200
R基础篇
矩阵(matrix)是一种特殊的向量,包含两个附加的属性:行数和列数。矩阵也和向量一样具有模式的概念。反过来向量不能看成只有一行或一列的矩阵,因为向量没有行数和列数这两个属性
数组(array)是R里更一般的对象,矩阵是数组的一个特殊形式。数组可以是多维的,例如一个三维数组可以包含行、列和层(layer),而一个矩阵只有行和列两个维度
数组或矩阵的内部存储都是按向量来存储,占用连续的存储空间
本文的写作顺序同样是讲述创建矩阵、矩阵的值操作、矩阵的运算以及一些方便的函数,最后介绍高维数组的创建和应用。数组可以看成是矩阵的推广,理解了矩阵就能更好的理解数组的概念了
> x <- matrix(c(1,2,3,4,5,6,7,8,9,10,11,12),nrow = 3,ncol = 4) #创建矩阵x,指定数据,行,列
> x
[,1] [,2] [,3] [,4]
[1,] 1 4 7 10
[2,] 2 5 8 11
[3,] 3 6 9 12
> (x <- matrix(c(1,2,3,4,5,6,7,8,9,10,11,12),nrow = 3)) #创建矩阵,指定数据、行数
[,1] [,2] [,3] [,4]
[1,] 1 4 7 10
[2,] 2 5 8 11
[3,] 3 6 9 12
> (x <- matrix(c(1,2,3,4,5,6,7,8,9,10,11,12),ncol = 4)) #创建矩阵,指定数据、列数
[,1] [,2] [,3] [,4]
[1,] 1 4 7 10
[2,] 2 5 8 11
[3,] 3 6 9 12
> (x <- matrix(c(1,2,3,4,5,6,7,8,9,10,11,12),nrow = 3,byrow = T)) #创建矩阵,按行排列
[,1] [,2] [,3] [,4]
[1,] 1 2 3 4
[2,] 5 6 7 8
[3,] 9 10 11 12
> (y <- matrix(nrow = 2, ncol = 2)) #创建一个2*2的空矩阵
[,1] [,2]
[1,] NA NA
[2,] NA NA
> y[1,1] <- 1 # 单独为矩阵中每个元素赋值
> y[1,2] <- 2
> y[2,1] <- 3
> y[2,2] <- 4
> y
[,1] [,2]
[1,] 1 2
[2,] 3 4
>
第二种是创建一个只有行和列的空矩阵,创建完成后再给矩阵中的每个元素赋值,这个时候就需要同时指定行数和列数;如最后一个例子所示
由于第二种方式比较麻烦,所以我们一般用第一种方式创建矩阵;
在设置数据向量的时候需要注意元素默认是按列来排列的,如果需要按行来排列,在matrix()函数后面加上byrow = T
参数;
使用byrow
参数时需要知道矩阵本身在内存中依然是按列存储的,参数byrow
只是改变了数据输入的顺序,当读取的数据文件是按行排列时,使用这个参数可能更方便些
取值
> (x <- matrix(c(1,2,3,4,5,6,7,8,9,10,11,12),nrow = 3,ncol = 4))
[,1] [,2] [,3] [,4]
[1,] 1 4 7 10
[2,] 2 5 8 11
[3,] 3 6 9 12
> x[2,3] # 取矩阵中第二行第三列元素的值
[1] 8
> x[2,] # 取矩阵的第二行
[1] 2 5 8 11
> x[,3] # 取矩阵的第三列
[1] 7 8 9
> class(x) # 矩阵x的类型是"matrix"
[1] "matrix"
> class(x[2,3]) # 矩阵中的每个元素的类型是一元向量
[1] "numeric"
> class(x[2,]) # 矩阵中取某一行结果得到向量
[1] "numeric"
> x[2] #取索引为2的元素
[1] 2
> x[4] #取索引为4的元素
[1] 4
> x[7] #取索引为7的元素
[1] 7
矩阵取值也是用[]
,括号其中的参数表示索引,取某个元素值时需要指定该元素的行号索引和列号索引,这俩属性用,
隔开,如x[2,3]
;取某行或者某列时,只需要指定该行或者该列就行,另一属性可以空着,但是中间的,
不能省,表示这个属性中的所有索引都取,如 x[2,]
和x[,3]
用class()函数查看对象的类型,可看成x是一个矩阵,按照线性代数中矩阵操作的逻辑,取矩阵的某一行或某一列得到的应该是1*n或者m*1的行矩阵,但是上述例子结果显示得到的却是一个向量,这就是R中会自动降维处理
从最后三个例子可以看出,矩阵本质上还是向量,并按列排列,当[]中只有一个索引值时,按向量的方式去取值
矩阵索引与修改值
> (x <- matrix(c(1,2,3,4,5,6,7,8,9,10,11,12),nrow = 3,ncol = 4))
[,1] [,2] [,3] [,4]
[1,] 1 4 7 10
[2,] 2 5 8 11
[3,] 3 6 9 12
> x[,2:3] #创建第二列到第三列的索引并按索引取值
[,1] [,2]
[1,] 4 7
[2,] 5 8
[3,] 6 9
> x[2:3,] #创建第二行到第三行的索引并按索引取值
[,1] [,2] [,3] [,4]
[1,] 2 5 8 11
[2,] 3 6 9 12
> x[2:3,2] #创建第二行到第三行,第二列的索引并按索引取值
[1] 5 6
> x[-1,] #去掉第一行后显示
[,1] [,2] [,3] [,4]
[1,] 2 5 8 11
[2,] 3 6 9 12
> x[,-1:-2] #去掉第一列到第二列后显示
[,1] [,2]
[1,] 7 10
[2,] 8 11
[3,] 9 12
> x[-1,-3] #同时去掉第一行与第三列的元素
[,1] [,2] [,3]
[1,] 2 5 11
[2,] 3 6 12
> x #取值只是取的原对象的副本,对原对象没有任何改变
[,1] [,2] [,3] [,4]
[1,] 1 4 7 10
[2,] 2 5 8 11
[3,] 3 6 9 12
> x[c(1,3),] <- matrix(c(1,1,8,12),nrow = 2) #生成2*2的矩阵,并按对应索引位置修改值
> x
[,1] [,2] [,3] [,4]
[1,] 1 8 1 8
[2,] 2 5 8 11
[3,] 1 12 1 12
> x[2:3,2:3] <- matrix(c(-1,-2,15,15),nrow = 2)
> x
[,1] [,2] [,3] [,4]
[1,] 1 8 1 8
[2,] 2 -1 15 11
[3,] 1 -2 15 12
> x[1:2,1:2] <- 0 #一元向量循环补齐后给矩阵对应索引的元素赋值
> x
[,1] [,2] [,3] [,4]
[1,] 0 0 1 8
[2,] 0 0 15 11
[3,] 1 -2 15 12
>
在x[1:2,1:2] <- 0
这个例子中,1:2,1:2
生成2*2的索引矩阵,然后一元向量0
会自动循环补齐成2*2的元素值全为0的矩阵,最后这个零矩阵给矩阵x对应索引处赋值从而得到结果
矩阵同向量一样,当索引为负值时,表示剔除掉该索引的行或列显示(这里说的剔除并不是删除对象中的元素的意思,只是显示的时候不显示剔除的元素);注意,如果指定的行索引和列索引都为负值,并不是单单剔除那一个元素,而是会同时剔除该行索引和列索引所在的整行和整列,如x[-1,-3]
矩阵元素筛选
> (x <- matrix(c(1,2,3,4,5,6,7,8,9,10,11,12),nrow = 3,ncol = 4))
[,1] [,2] [,3] [,4]
[1,] 1 4 7 10
[2,] 2 5 8 11
[3,] 3 6 9 12
> x[,2] >= 5
[1] FALSE TRUE TRUE
> x[x[,2] >= 5] # 结果展开成列向量
[1] 2 3 5 6 8 9 11 12
> x[x[,2] >= 5, ] #结果是矩阵
[,1] [,2] [,3] [,4]
[1,] 2 5 8 11
[2,] 3 6 9 12
x[,2] >= 5
会由一元向量5
循环补齐成5 5 5
,x[,2]
为向量4 5 6
与5 5 5
对应元素做>=
运算,得到索引向量FALSE TRUE TRUE
x[x[,2] >= 5]
表示矩阵x取第二行和第三行索引对应的元素,并按列排列生成向量
而x[x[,2] >= 5, ]
表示矩阵x取第二行和第三行索引对应的元素,并生成矩阵
增加或删除矩阵的行或列
> (x <- matrix(c(1,2,3,4,5,6,7,8,9,10,11,12),nrow = 3,ncol = 4))
[,1] [,2] [,3] [,4]
[1,] 1 4 7 10
[2,] 2 5 8 11
[3,] 3 6 9 12
> (y <- c(-1,-2,-3))
[1] -1 -2 -3
> rbind(y,x) #向量y与矩阵x按行组合
[,1] [,2] [,3] [,4]
y -1 -2 -3 -1
1 4 7 10
2 5 8 11
3 6 9 12
Warning message:
In rbind(y, x) :
number of columns of result is not a multiple of vector length (arg 1)
> cbind(y,x) #向量y与矩阵x按列组合
y
[1,] -1 1 4 7 10
[2,] -2 2 5 8 11
[3,] -3 3 6 9 12
> (x <- x[-2,]) #删除一行
[,1] [,2] [,3] [,4]
[1,] 1 4 7 10
[2,] 3 6 9 12
>
改变矩阵的大小可以通过rbind()函数和cbind()函数来实现,rbind()是将向量或者矩阵按行组合生成新的矩阵,cbind()是将向量或矩阵按列组合生成新的矩阵,如果两个对象的行数或列数不等时,R会自动循环补齐
以rbind(y,x)
为例,向量y与矩阵x按行绑定,但是y只有3列,x有4列,于是y先循环补齐成-1 -2 -3 -1
,又因为参数y在前面,所以在与x按行组合时在矩阵x的上方,从而得到结果矩阵
事实上,矩阵的长度和维度在矩阵创建时就已经固定,因此不能增加或删除行或列,所以矩阵的增加或删除同向量一样,并不是在原来的矩阵上增加或删除,而是重新生成了一个矩阵,然后让变量绑定到新矩阵上,给人的感觉就好像是在原来的矩阵上增加或删除行或列
对矩阵行和列调用函数
> (x <- matrix(c(1,2,3,4,5,6,7,8,9,10,11,12),nrow = 3,ncol = 4))
[,1] [,2] [,3] [,4]
[1,] 1 4 7 10
[2,] 2 5 8 11
[3,] 3 6 9 12
> apply(x,2,mean) #对x按每列使用mean()函数
[1] 2 5 8 11
> apply(x,1,mean) #对x按每行使用mean()函数
[1] 5.5 6.5 7.5
> multiple <- function(x) x*2 #创建个函数,作用是将x加倍,可暂时不用管函数如何创建
> (y <- apply(x,1,multiple)) #将x每一行加倍
[,1] [,2] [,3]
[1,] 2 4 6
[2,] 8 10 12
[3,] 14 16 18
[4,] 20 22 24
> t(y) # t()函数可将矩阵转置
[,1] [,2] [,3] [,4]
[1,] 2 8 14 20
[2,] 4 10 16 22
[3,] 6 12 18 24
>
apply(m,dimcode,f,fargs)函数系列是R中最受欢迎也是最常用的函数,该函数系列包括apply()、tapply()和lapply()。
该函数的参数列表中m是一个矩阵;dimcode是维度编号,若取值为1代表对每一行应用f函数,若取值为2代表对每一列应用f函数;f是一个函数,应用在行或列上;fargs是f的可选参数集;apply()函数所得的结果都会按照列排列,如果想要得到按行排列的结果矩阵,可以使用 t() 函数进行转置
上例中以apply(x,2,mean)
来说明,mean()函数是求向量的均值,按照apply函数的规则,对矩阵x的每一列应用mean()函数,则第一列列表项1 2 3
应用mean()函数求得均值为 (1+2+3)/3=2,第二列4 5 6
求均值得(4+5+6)/3=5,第三列7 8 9
求均值得(7+8+9)/3=8,第四列10 11 12
求均值得(10+11+12)/3=11,所以该表达式的结果为2 5 8 11
而在apply(x,1,multiple)
例子中,对x的每一行应用multiple()函数,作用是将x的每一行乘以2,即将每一行的每个元素加倍;可以看到得到的结果是4*3的矩阵,而不是3*4,这就是因为apply()函数的结果都是按列排列的
其中关于函数创建这个表达式可以暂时不用管,在我们学了自定义函数后应用于第三个参数,会再来使用apply()函数
向量与矩阵的差异
> z <- matrix(1:8,nrow=4)
> z
[,1] [,2]
[1,] 1 5
[2,] 2 6
[3,] 3 7
[4,] 4 8
> length(z) #矩阵也是一个向量,可以查看向量的长度
[1] 8
> class(z) #作为矩阵,不仅仅是一个向量,还具有矩阵类属性
[1] "matrix"
> attributes(z) #矩阵具有尺寸属性,行和列
$dim
[1] 4 2
> dim(z) #dim()函数可以查看矩阵的大小
[1] 4 2
> nrow(z) #nrow()函数可以查看矩阵的行数
[1] 4
> ncol(z) #ncol()函数可以查看矩阵的列数
[1] 2
>
矩阵是一个向量,只是多了两个属性:行数和列数。从面向对象编程的角度说,矩阵类(matrix class)是实际存在的。
避免意外降维
在统计学领域,“降维”(dimension reduction)是有益的,也存在很多降维的统计学方法,因为这样可以降低大量的计算。
但是在R中,降维指的完全是另外一回事,而且通常要避免,例如上面提取矩阵的某一行时:
> (x <- matrix(c(1,2,3,4,5,6,7,8,9,10,11,12),nrow = 3,ncol = 4))
[,1] [,2] [,3] [,4]
[1,] 1 4 7 10
[2,] 2 5 8 11
[3,] 3 6 9 12
> (y <- x[2,]) #取矩阵的第二行
[1] 2 5 8 11
> attributes(x) # 查看矩阵x的属性
$dim
[1] 3 4
> attributes(y) # 查看子集y的属性为NULL
NULL
> str(x) # 打印x对象可看成是矩阵
num [1:3, 1:4] 1 2 3 4 5 6 7 8 9 10 ...
> str(y) # 打印y对象可看成是一个向量
num [1:4] 2 5 8 11
> (z <- x[2, , drop = FALSE]) #使用drop参数可以防止意外降维
[,1] [,2] [,3] [,4]
[1,] 2 5 8 11
> dim(z) #此时查看子集z对象具有尺寸属性
[1] 1 4
>
上述例子显示,当直接提取矩阵的某一行时,结果看着没问题,但是结果对象的类型就不对了,一般情况下我们想要的是从矩阵中提取一行,得到的结果应该是1*n的矩阵,而不是向量,这就是R中的意外降维;
一般情况下将一行或一列子集看成是向量也没问题,但是在有些大量矩阵操作的程序中会出错;可以使用drop=FALSE
参数来防止意外降维
矩阵的行列命名
> (x <- matrix(c(1,2,3,4,5,6,7,8,9,10,11,12),nrow = 3,ncol = 4))
[,1] [,2] [,3] [,4]
[1,] 1 4 7 10
[2,] 2 5 8 11
[3,] 3 6 9 12
> colnames(x) #查看矩阵的列名
NULL
> colnames(x) <- c("a","b") #为矩阵设置列名时,名称的长度需要与矩阵的列数一致
Error in `colnames<-`(`*tmp*`, value = c("a", "b")) :
'dimnames'的长度[2]必需与陈列范围相等
> colnames(x) <- c("a","b","c","d")
> x #可以看到矩阵的列名设置好了
a b c d
[1,] 1 4 7 10
[2,] 2 5 8 11
[3,] 3 6 9 12
> rownames(x) # 查看矩阵的行名称
NULL
> rownames(x) <- c("f","g","h") #为矩阵设置行名称
> x
a b c d
f 1 4 7 10
g 2 5 8 11
h 3 6 9 12
> x["g","c"] # 可通过行和列的名称来访问元素
[1] 8
>
colnames(x)函数可以给矩阵的列命名,也可以查看矩阵的列名
rownames(x)函数可以给矩阵的行命名,也可以查看矩阵的行名
访问矩阵元素最直接的方法是通过行号和列号,也可以使用行名和列名,但是一般在编写R代码时,给行和列命名并不那么重要,只是在分析某些数据时会有用,可以方便查看某行或某列的数据代表什么意思
> x <- matrix(1:12,nrow = 4)
> x
[,1] [,2] [,3]
[1,] 1 5 9
[2,] 2 6 10
[3,] 3 7 11
[4,] 4 8 12
> t(x)
[,1] [,2] [,3] [,4]
[1,] 1 2 3 4
[2,] 5 6 7 8
[3,] 9 10 11 12
t(x)函数可以对矩阵x做转置变换
> rowSums(x) # 对矩阵x按行求每行的和
[1] 15 18 21 24
> rowMeans(x) # 对矩阵x按行求每行的均值
[1] 5 6 7 8
> colSums(x) #对矩阵x按列求每列的和
[1] 10 26 42
> colMeans(x) #对矩阵x按列求每列的均值
[1] 2.5 6.5 10.5
这四个函数分别于apply(x,1,sum),apply(x,1,mean),apply(x,2,sum),apply(x,2,mean)效果一样
> y <- matrix(rnorm(9),nrow = 3,ncol = 3)
> y
[,1] [,2] [,3]
[1,] 0.3683004 0.06304847 0.9411764
[2,] -0.3673749 0.53875845 1.0974140
[3,] 2.0538858 0.99482935 1.7670388
> det(y)
[1] -1.25386
det(y)可以求方阵y的行列式,矩阵y必须是行列相等的方阵
> (A = B = matrix(1:12,nrow = 3,ncol = 4))
[,1] [,2] [,3] [,4]
[1,] 1 4 7 10
[2,] 2 5 8 11
[3,] 3 6 9 12
> A+B #矩阵加法是按对应索引元素相加
[,1] [,2] [,3] [,4]
[1,] 2 8 14 20
[2,] 4 10 16 22
[3,] 6 12 18 24
> A - B #矩阵减法是按对应索引元素相减
[,1] [,2] [,3] [,4]
[1,] 0 0 0 0
[2,] 0 0 0 0
[3,] 0 0 0 0
> c <- 2
> c*A #矩阵数乘是对矩阵中的每个元素乘以这个数
[,1] [,2] [,3] [,4]
[1,] 2 8 14 20
[2,] 4 10 16 22
[3,] 6 12 18 24
矩阵元素的加减是指维数相同的矩阵,处于同行和同列位置的元素进行加减;矩阵数乘是指一个常数与一个矩阵相乘,其结果是矩阵中的每个元素都与这个常数相乘
> (A <- matrix(1:12,nrow = 3,ncol = 4))
[,1] [,2] [,3] [,4]
[1,] 1 4 7 10
[2,] 2 5 8 11
[3,] 3 6 9 12
> (B <- matrix(1:12,nrow = 4,ncol = 3))
[,1] [,2] [,3]
[1,] 1 5 9
[2,] 2 6 10
[3,] 3 7 11
[4,] 4 8 12
> A %*% B #矩阵乘法
[,1] [,2] [,3]
[1,] 70 158 246
[2,] 80 184 288
[3,] 90 210 330
A是一个3*4的矩阵,B是一个4*3的矩阵,则A与B要做矩阵乘法,需要使用%*%
符号,可以生成一个3*3的矩阵;m*k的矩阵与k*n的矩阵做矩阵乘法可以得到一个m*n的矩阵,由于矩阵乘法不满足交换律,所以k*n的矩阵与m*k的矩阵相乘没有意义
如果要计算A'B,A为n*m的矩阵,B为n*k的矩阵,A'表示矩阵A的转置,则可以用t()函数配合 %*% 来实现
> (A = B = matrix(1:12,nrow = 4))
[,1] [,2] [,3]
[1,] 1 5 9
[2,] 2 6 10
[3,] 3 7 11
[4,] 4 8 12
> t(A) %*% B
[,1] [,2] [,3]
[1,] 30 70 110
[2,] 70 174 278
[3,] 110 278 446
> crossprod(A,B)
[,1] [,2] [,3]
[1,] 30 70 110
[2,] 70 174 278
[3,] 110 278 446
>
先对矩阵A转置,然后求转置后的矩阵与B做矩阵乘法,也可以使用crossprod(A,B)函数来实现
> (A <- matrix(1:12, 4, 3))
[,1] [,2] [,3]
[1,] 1 5 9
[2,] 2 6 10
[3,] 3 7 11
[4,] 4 8 12
> (C <- matrix(1:4,2,2))
[,1] [,2]
[1,] 1 3
[2,] 2 4
> kronecker(A,C)
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1 3 5 15 9 27
[2,] 2 4 10 20 18 36
[3,] 2 6 6 18 10 30
[4,] 4 8 12 24 20 40
[5,] 3 9 7 21 11 33
[6,] 6 12 14 28 22 44
[7,] 4 12 8 24 12 36
[8,] 8 16 16 32 24 48
>
矩阵的Kronecker积可以用函数kronecher(A,C)来计算,n*m的矩阵与h*k的矩阵做Kronecker积,结果是一个nk*mk维矩阵;矩阵A中的每个元素与矩阵B做数乘,然后按照这个元素在矩阵A中的位置排列生成新矩阵
矩阵运算的数学原理现在不熟悉也没关系,在基础篇结束后我们会用专门的篇幅介绍
> (A <- matrix(rnorm(9),nrow = 3,ncol = 3))
[,1] [,2] [,3]
[1,] -1.160287 -2.5451521 -0.1679758
[2,] 0.702371 1.5207645 -1.0167863
[3,] -2.004842 0.5957342 -0.4857430
> solve(A)
[,1] [,2] [,3]
[1,] 0.02050442 0.20607670 -0.438462771
[2,] -0.36696295 -0.03497988 0.200122201
[3,] -0.53468711 -0.89345607 -0.003564942
> A %*% solve(A)
[,1] [,2] [,3]
[1,] 1 -2.775558e-17 4.857226e-17
[2,] 0 1.000000e+00 -3.165870e-17
[3,] 0 0.000000e+00 1.000000e+00
> round(A %*% solve(A)) #使用round()函数后可以更好的得到结果
[,1] [,2] [,3]
[1,] 1 0 0
[2,] 0 1 0
[3,] 0 0 1
>
solve(A,B)函数可以求解方程组,AX=B;当第二个参数缺省时,可以求解矩阵A的逆矩阵
我们知道一个矩阵与该矩阵的逆矩阵做矩阵乘法,会得到一个单位矩阵,从而可以验证solve(A)函数可以求矩阵A的逆矩阵
> (A <- matrix(1:9, nrow = 3))
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9
> Aeigen <- eigen(A)
> Aeigen
$values
[1] 1.611684e+01 -1.116844e+00 -5.700691e-16
$vectors
[,1] [,2] [,3]
[1,] -0.4645473 -0.8829060 0.4082483
[2,] -0.5707955 -0.2395204 -0.8164966
[3,] -0.6770438 0.4038651 0.4082483
>
在线性代数中一个非常重要的矩阵操作是求矩阵的特征值和特征向量,在实际应用中求特征值有着非常广泛的用途
可以通过对矩阵A进行谱分解来得到矩阵的特征值和特征向量。矩阵A的谱分解为:A = UVU',其中U的列为A的特征值所对应的特征向量
在R中可以用eigen()函数得到U和V,结果中Aeigen$values
即为特征值,Aeigen$vectors
即为特征向量
> (A <- matrix(1:18, 3, 6))
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1 4 7 10 13 16
[2,] 2 5 8 11 14 17
[3,] 3 6 9 12 15 18
> (y <- svd(A))
$d
[1] 4.589453e+01 1.640705e+00 1.366522e-15
$u
[,1] [,2] [,3]
[1,] -0.5290354 0.74394551 0.4082483
[2,] -0.5760715 0.03840487 -0.8164966
[3,] -0.6231077 -0.66713577 0.4082483
$v
[,1] [,2] [,3]
[1,] -0.07736219 -0.71960032 -0.4076688
[2,] -0.19033085 -0.50893247 0.5745647
[3,] -0.30329950 -0.29826463 -0.0280114
[4,] -0.41626816 -0.08759679 0.2226621
[5,] -0.52923682 0.12307105 -0.6212052
[6,] -0.64220548 0.33373889 0.2596585
> y$u %*% diag(y$d) %*% t(y$v)
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1 4 7 10 13 16
[2,] 2 5 8 11 14 17
[3,] 3 6 9 12 15 18
>
线性代数中有两个经常用到的矩阵分解,其中之一就是奇异值分解(SVD分解)
A为m*n矩阵,矩阵的秩为r,A可以分解为A=UDV',其中U'U = V'V = 1。
在R中可以用svd(A)函数求U、D、V;结果可以通过UDV' = A 来检验
> (A <- matrix(1:12, 4, 3))
[,1] [,2] [,3]
[1,] 1 5 9
[2,] 2 6 10
[3,] 3 7 11
[4,] 4 8 12
> qr(A) # 对矩阵A进行QR分解
$qr
[,1] [,2] [,3]
[1,] -5.4772256 -12.7801930 -2.008316e+01
[2,] 0.3651484 -3.2659863 -6.531973e+00
[3,] 0.5477226 -0.3781696 1.601186e-15
[4,] 0.7302967 -0.9124744 -5.547002e-01
$rank
[1] 2
$qraux
[1] 1.182574 1.156135 1.832050
$pivot
[1] 1 2 3
attr(,"class")
[1] "qr"
>
线性代数中常用的另外一种矩阵分解就是QR分解;A为m*n矩阵可以进行QR分解:A = QR,其中Q'Q = 1;在R中可以用函数qr(A)来完成分解,结果中$rank
项返回的是矩阵的秩,$qr
项包含了Q矩阵和R矩阵的信息,可以分别通过qr.Q()函数和qr.R()函数查看:
> qr.Q(qr(A))
[,1] [,2] [,3]
[1,] -0.1825742 -8.164966e-01 -0.4000874
[2,] -0.3651484 -4.082483e-01 0.2546329
[3,] -0.5477226 -1.665335e-16 0.6909965
[4,] -0.7302967 4.082483e-01 -0.5455419
> qr.R(qr(A))
[,1] [,2] [,3]
[1,] -5.477226 -12.780193 -2.008316e+01
[2,] 0.000000 -3.265986 -6.531973e+00
[3,] 0.000000 0.000000 1.601186e-15
>
这里只列举出了一些常用的矩阵运算,还有其他一些高级的运算就不一一细讲,感兴趣的读者可以借助网络自行学习
矩阵是二维的数据结构,能表示两个维度的数据,但是当我们想表示更高维度的数据时,就需要用到数组(arrays)了
举个简单的例子,考虑学生和考试成绩的数据。假设每次考试分两部分,因此每次考试需要给每个学生记录两个分数。假设有两次考试,只有三个学生。其中第一次考试的数据为:
> (firsttest <- matrix(c(46,30,21,25,50,50),nrow=3,byrow=T))
[,1] [,2]
[1,] 46 30
[2,] 21 25
[3,] 50 50
第二次考试的成绩是:
> (secondtest <- matrix(c(46,43,41,35,50,50),nrow=3,byrow=T))
[,1] [,2]
[1,] 46 43
[2,] 41 35
[3,] 50 50
现在想把两次考试的成绩合并到一个数据结构里,命名为tests。则tests应该分两个数据层,一层对应一次考试,每层有三行两列。若假设firsttest在第一层,则用array()函数创建这样的数据结构为:
> (tests <- array(data = c(firsttest,secondtest),dim = c(3,2,2)))
, , 1
[,1] [,2]
[1,] 46 30
[2,] 21 25
[3,] 50 50
, , 2
[,1] [,2]
[1,] 46 43
[2,] 41 35
[3,] 50 50
>
其中参数dim = c(3,2,2)
指明数据共有两层(第三个参数),每层分别有三行(第一个参数)两列(第二个参数),这个参数最后会成为数组的dim属性:
> attributes(tests)
$dim
[1] 3 2 2
array数组的每个元素现在都有三个下标,比矩阵多一个,访问数组中的元素时需要按照dim属性的顺序设置索引:
> tests[1,2,1] #查看第一个同学第一次考试的第二门成绩
[1] 30
>
我们可以把两个矩阵合并成一个三维数组,同样也可以把两个或多个三维数组合并成四维数组,以此类推