2.5 多维数组和矩阵

2.5.1 生成数组或矩阵

数组有一个特征属性叫做维数向量(dim属性),维数向量是一个元素取正整数的向量,其长度是数组的维数,比如维数向量有两个元素时数组为2维数组(矩阵)。维数向量的每一个元素指定了该下标的上界,下标的下界总为1

1.将向量定义成数组

向量只有定义了维数向量(dim属性)后才能被看作是数组
> z<-1:12
> dim(z)<-c(3,4);z
     [,1] [,2] [,3] [,4]
[1,]    1    4    7   10
[2,]    2    5    8   11
[3,]    3    6    9   12
> dim(z)<-12;z
 [1]  1  2  3  4  5  6  7  8  9 10 11 12

2.用array()函数构造多维数组

R软件可以用array()函数直接构造数组,其构造形式为
array(data=NA,dim=length(data),dimnames=NULL)
其中data是一个向量数据,dim是数组各维的长度,缺省时是原向量的长度,dimnames是数组维的名字,缺省时为空
> X<- array(1:20,dim=c(4,5));X
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    5    9   13   17
[2,]    2    6   10   14   18
[3,]    3    7   11   15   19
[4,]    4    8   12   16   20
> Z<- array(0,dim=c(3,4,2));Z
, , 1

     [,1] [,2] [,3] [,4]
[1,]    0    0    0    0
[2,]    0    0    0    0
[3,]    0    0    0    0

, , 2

     [,1] [,2] [,3] [,4]
[1,]    0    0    0    0
[2,]    0    0    0    0
[3,]    0    0    0    0

3.用matrix()函数构造矩阵

函数matrix()是构造矩阵的函数,其构造形式为
matrix(data=NA,nrow=1,ncol=1,byrow=FALSE,dimnames=NULL)
其中data是一个向量数据,nrow是矩阵的行数,ncol是矩阵的列数。当byrow=TRUE时,生成矩阵的数据按行放置,缺省时相当于byrow=FALSE,数据按列放置。dimnames是数组维数的名字,缺省时为空
> A<- matrix(1:15,nrow=3,ncol=5,byrow = TRUE);A
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    2    3    4    5
[2,]    6    7    8    9   10
[3,]   11   12   13   14   15
> A<- matrix(1:15,nrow=3,byrow=TRUE);A
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    2    3    4    5
[2,]    6    7    8    9   10
[3,]   11   12   13   14   15

2.5.2 数组下标

1.数组下标

> a<-1:24
> dim(a)<-c(2,3,4)
> a[2,1,2]
[1] 8
> a[1,2:3,2:3]
     [,1] [,2]
[1,]    9   15
[2,]   11   17
> a
, , 1

     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6

, , 2

     [,1] [,2] [,3]
[1,]    7    9   11
[2,]    8   10   12

, , 3

     [,1] [,2] [,3]
[1,]   13   15   17
[2,]   14   16   18

, , 4

     [,1] [,2] [,3]
[1,]   19   21   23
[2,]   20   22   24

> a[1,,]
     [,1] [,2] [,3] [,4]
[1,]    1    7   13   19
[2,]    3    9   15   21
[3,]    5   11   17   23
> a[,2,]
     [,1] [,2] [,3] [,4]
[1,]    3    9   15   21
[2,]    4   10   16   22
> a[1,1,]
[1]  1  7 13 19
> a[3:10]
[1]  3  4  5  6  7  8  9 10

2.不规则的数组下标

在R语言中,甚至可以把数组中的任意位置的元素作为数组访问,其方法是用一个二维数组作为数组的下标,二维数组的每一个行是一个元素的下标,列数为数组的维数。
例如,要把上面的形状为234的数组a的第[1,1,1],[2,2,3],[1,3,4],[2,1,4]号共四个元素作为一个整体访问,先定义一个包含这些下标作为行的二维数组:
> b<- matrix(c(1,1,1,2,2,3,1,3,4,2,1,4),ncol=3,byrow = TRUE);b
     [,1] [,2] [,3]
[1,]    1    1    1
[2,]    2    2    3
[3,]    1    3    4
[4,]    2    1    4
> a[b]
[1]  1 16 23 20
> a[b]<-c(101,102,103,104);a[b]
[1] 101 102 103 104
> a
, , 1

     [,1] [,2] [,3]
[1,]  101    3    5
[2,]    2    4    6

, , 2

     [,1] [,2] [,3]
[1,]    7    9   11
[2,]    8   10   12

, , 3

     [,1] [,2] [,3]
[1,]   13   15   17
[2,]   14  102   18

, , 4

     [,1] [,2] [,3]
[1,]   19   21  103
[2,]  104   22   24

2.5.3 数组的四则运算

> A<- matrix(1:6,nrow = 2,byrow = TRUE);A
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    4    5    6
> B<- matrix(1:6,nrow = 2);B
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6
> C<- matrix(c(1,2,2,3,3,4),nrow = 2);C
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    2    3    4
> D<- 2*C+A/B;D
     [,1]     [,2] [,3]
[1,]    3 4.666667  6.6
[2,]    6 7.250000  9.0
形状不一致的向量(或数组)也可以进行四则运算,一般的规则是将向量(或数组)中的数据与对应向量(或数组)中的数据进行运算,把短向量(或数组)的数据循环使用,从而可以与长向量(或数组)数据进行匹配,并尽可能保留共同的数组属性
> x1<-c(100,200)
> x2<-1:6
> x1+x2
[1] 101 202 103 204 105 206
> x3<- matrix(1:6,nrow = 3)
> x1+x3
     [,1] [,2]
[1,]  101  204
[2,]  202  105
[3,]  103  206
> x2<- 1:5
> x1+x2
[1] 101 202 103 204 105
Warning message:
In x1 + x2 : 长的对象长度不是短的对象长度的整倍数

2.5.4 矩阵的运算

1.转置运算

> A<-matrix(1:6,nrow=2);A
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6
> t(A)
     [,1] [,2]
[1,]    1    2
[2,]    3    4
[3,]    5    6
> ### 2.求方阵的行列式
> det(matrix(1:4,ncol=2))
[1] -2
> ### 3.向量的内积
> x<-1:5;y<-2*1:5
> x%*%y
     [,1]
[1,]  110
函数crossprod()是内积运算函数(表示交叉乘积),crossprod(x,y)计算向量x与y的内积,即't(x) %*% y'
> crossprod(x,y)
     [,1]
[1,]  110
> x;t(x)
[1] 1 2 3 4 5
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    2    3    4    5
crossprod(x)表示x与x的内积,即\(||x||^2_2\)
类似地,tcrossprod(x,y)表示'x%*% t(y)',即x与y的外积,也称为叉积。tcrossprod(x)表示x与x作外积
> tcrossprod(x,y)
     [,1] [,2] [,3] [,4] [,5]
[1,]    2    4    6    8   10
[2,]    4    8   12   16   20
[3,]    6   12   18   24   30
[4,]    8   16   24   32   40
[5,]   10   20   30   40   50
> tcrossprod(x)
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    2    3    4    5
[2,]    2    4    6    8   10
[3,]    3    6    9   12   15
[4,]    4    8   12   16   20
[5,]    5   10   15   20   25

4.向量的外积(叉积)

设x,y是n维向量,则x %o% y表示x与y作外积
> x<-1:5;y<-2*1:5
> x %o% y
     [,1] [,2] [,3] [,4] [,5]
[1,]    2    4    6    8   10
[2,]    4    8   12   16   20
[3,]    6   12   18   24   30
[4,]    8   16   24   32   40
[5,]   10   20   30   40   50
outer()是外积运算函数,outer(x,y)计算向量x与y的外积,它等价于 x %o% y
函数outer()的一般调用格式为
outer(X,Y,fun='*',...)
其中X,Y矩阵(或向量),fun是作外积运算函数,缺省值为乘法运算,函数outer()在绘制三维曲面时非常有用,它可以生成一个X和Y的网格,
> outer(x,y)
     [,1] [,2] [,3] [,4] [,5]
[1,]    2    4    6    8   10
[2,]    4    8   12   16   20
[3,]    6   12   18   24   30
[4,]    8   16   24   32   40
[5,]   10   20   30   40   50

5.矩阵的乘法

如果矩阵A和B具有相同的维数,则AB表示矩阵中对应的元素的乘积,A %% B表示通常意义下的两个矩阵的乘积(要求矩阵A的列数等于矩阵B的行数)
> A<-array(1:9,dim=(c(3,3)))
> B<-array(9:1,dim=(c(3,3)))
> C<-A*B;C
     [,1] [,2] [,3]
[1,]    9   24   21
[2,]   16   25   16
[3,]   21   24    9
> D<- A%*%B;D
     [,1] [,2] [,3]
[1,]   90   54   18
[2,]  114   69   24
[3,]  138   84   30
crossprod(A,B)表示的是t(A)%%B,函数tcrossprod(A,B)表示的是A%%t(B)
> A
     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9
> B<-c(1:3);B
[1] 1 2 3
> crossprod(B,A)
     [,1] [,2] [,3]
[1,]   14   32   50
> t(B);crossprod(t(B),A)
     [,1] [,2] [,3]
[1,]    1    2    3
Error in crossprod(t(B), A) : 非整合参数
> tcrossprod(A,B)
Error in tcrossprod(A, B) : 非整合参数
> tcrossprod(A,t(B))
     [,1]
[1,]   30
[2,]   36
[3,]   42

6.生成对角阵和矩阵取对角运算

函数diag()依赖于它的变量,当v是一个向量时,diag(v)表示以v的元素为对角线元素的对角阵。当M是一个矩阵时,则diag(M)表示的是取M对角线上的元素的向量
> v<-c(1,4,5)
> diag(v)
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    4    0
[3,]    0    0    5
> M<-array(1:9,dim=c(3,3));M
     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9
> diag(M)
[1] 1 5 9

7.解线性方程组和求矩阵的逆矩阵

若求解线性方程组Ax=b,其命令形式为solve(A,b),求矩阵A的逆,其命令形式为solve(A)
> A<-t(array(c(1:8,10),dim=c(3,3)))
> b<-c(1,1,1)
> x<-solve(A,b);x
[1] -1.000000e+00  1.000000e+00  3.330669e-16
> B<-solve(A);B
           [,1]      [,2] [,3]
[1,] -0.6666667 -1.333333    1
[2,] -0.6666667  3.666667   -2
[3,]  1.0000000 -2.000000    1
> A %*% B
     [,1]         [,2]          [,3]
[1,]    1 8.881784e-16 -4.440892e-16
[2,]    0 1.000000e+00 -1.776357e-15
[3,]    0 0.000000e+00  1.000000e+00

8.求矩阵的特征值与特征值向量

eigen(Sm)是求对称矩阵Sm的特征值与特征向量,其命令形式为
ev<-eigen(Sm)
ev存放着对称矩阵Sm特征值和特征向量,是由列表形式给出,其中ev\(values是Sm的特征值构成的向量,ev\)vectors是Sm的特征向量构成的矩阵
> A
     [,1] [,2] [,3]
[1,]    1    2    3
[2,]    4    5    6
[3,]    7    8   10
> Sm<-crossprod(A,A);Sm
     [,1] [,2] [,3]
[1,]   66   78   97
[2,]   78   93  116
[3,]   97  116  145
> ev<-eigen(Sm);ev
eigen() decomposition
$values
[1] 303.19533618   0.76590739   0.03875643

$vectors
           [,1]         [,2]       [,3]
[1,] -0.4646675  0.833286355  0.2995295
[2,] -0.5537546 -0.009499485 -0.8326258
[3,] -0.6909703 -0.552759994  0.4658502

9.矩阵的奇异值分解

函数svd(A)是对称A作奇异值分解,即A=UDV',其中U,V是正交阵,D为对角阵,也就是矩阵A的奇异值。svd(A)的返回值也是列表,svd(A)\(d表示矩阵A的奇异值,即矩阵D的对角线上的元素,svd(A)\)u对应的是正交阵U,svd(A)$v对应的是正交阵V
> svdA<-svd(A);svdA
$d
[1] 17.4125052  0.8751614  0.1968665

$u
           [,1]        [,2]       [,3]
[1,] -0.2093373  0.96438514  0.1616762
[2,] -0.5038485  0.03532145 -0.8630696
[3,] -0.8380421 -0.26213299  0.4785099

$v
           [,1]         [,2]       [,3]
[1,] -0.4646675 -0.833286355  0.2995295
[2,] -0.5537546  0.009499485 -0.8326258
[3,] -0.6909703  0.552759994  0.4658502

> attach(svdA)
The following object is masked _by_ .GlobalEnv:

    v

> u %*% diag(d) %*% t(v)
Error in u %*% diag(d) %*% t(v) : 非整合参数

10.求矩阵行列式的值

> det(A)
[1] -3

11.最小拟合与QR分解

> x<-c(0.0,0.2,0.4,0.6,0.8)
> y<-c(0.9,1.9,2.8,3.3,4.2)
> lsfit.sol<-lsfit(x,y);lsfit.sol
$coefficients
Intercept         X 
     1.02      4.00 

$residuals
[1] -0.12  0.08  0.18 -0.12 -0.02

$intercept
[1] TRUE

$qr
$qt
[1] -5.85849810  2.52982213  0.23749843 -0.02946714  0.10356728

$qr
      Intercept          X
[1,] -2.2360680 -0.8944272
[2,]  0.4472136  0.6324555
[3,]  0.4472136 -0.1954395
[4,]  0.4472136 -0.5116673
[5,]  0.4472136 -0.8278950

$qraux
[1] 1.447214 1.120788

$rank
[1] 2

$pivot
[1] 1 2

$tol
[1] 1e-07

attr(,"class")
[1] "qr"
\(coefficients是拟合系数,\)residuals是拟合残差
与lsfit()函数有密切关系的函数是ls.diag(),它给出拟合的进一步统计信息
另一个最小二乘拟合有密切关系的函数是QR分解函数qr(),qe.coef(),qr.fitted()和qr.resid()
> X<-matrix(c(rep(1,5),x),ncol = 2);X
     [,1] [,2]
[1,]    1  0.0
[2,]    1  0.2
[3,]    1  0.4
[4,]    1  0.6
[5,]    1  0.8
> Xplus<-qr(X);Xplus
$qr
           [,1]       [,2]
[1,] -2.2360680 -0.8944272
[2,]  0.4472136  0.6324555
[3,]  0.4472136 -0.1954395
[4,]  0.4472136 -0.5116673
[5,]  0.4472136 -0.8278950

$rank
[1] 2

$qraux
[1] 1.447214 1.120788

$pivot
[1] 1 2

attr(,"class")
[1] "qr"
QR分解函数qr()输入的设计矩阵需要加以1为元素的列,其返回值是列表,其中\(qr矩阵的上三角阵是QR分解中得到的R矩阵,下三角阵是QR分解得到的正交阵Q的部分信息,\)qraux是Q的附加信息
可用QR分解得到的结果计算最小二乘的系数
> b<- qr.coef(Xplus,y);b
[1] 1.02 4.00
得到的系数与函数lsfit()得到的结果相同,但是为什么用这种方法计算呢?这是因为用QR分解在计算最小二乘拟合时,其计算误差比一般方法要小
类似的,可以用QR分解得到最小二乘的拟合值和残差值
> fit<-qr.fitted(Xplus,y);fit
[1] 1.02 1.82 2.62 3.42 4.22
> res<-qr.resid(Xplus,y);res
[1] -0.12  0.08  0.18 -0.12 -0.02

2.5.5 与矩阵(数组)运算相关的函数

1.取矩阵的维数

> A<-matrix(1:6,nrow=2);A
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6
> dim(A)
[1] 2 3
> nrow(A)
[1] 2
> ncol(A)
[1] 3

2.矩阵的合并

函数cbind()把其自变量横向拼成一个大矩阵,rbind()把其自变量纵向拼成一个大矩阵
> x1<-rbind(c(1,2),c(3,4));x1
     [,1] [,2]
[1,]    1    2
[2,]    3    4
> x2<-10+x1
> x3<-cbind(x1,x2);x3
     [,1] [,2] [,3] [,4]
[1,]    1    2   11   12
[2,]    3    4   13   14
> x4<-rbind(x1,x2);x4
     [,1] [,2]
[1,]    1    2
[2,]    3    4
[3,]   11   12
[4,]   13   14
> cbind(1,x1)
     [,1] [,2] [,3]
[1,]    1    1    2
[2,]    1    3    4

3.矩阵的拉直

> A<-matrix(1:6,nrow=2);A
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6
> as.vector(A)
[1] 1 2 3 4 5 6

4.数组的维名字

> X<-matrix(1:6,ncol=2,dimnames = list(c("one","two","three"),c("First","Second")),byrow=T);X
      First Second
one       1      2
two       3      4
three     5      6
> X<-matrix(1:6,ncol=2,byrow=T)
> dimnames(X)<-list(c("one","two","three"),c("First","Second"))
> colnames(X)
[1] "First"  "Second"
> rownames(X)
[1] "one"   "two"   "three"

5.矩阵的广义转置

可以用aperm(A,perm)函数把数组A的各维按perm中指定的新次序重新排列
> A<-array(1:24,dim=c(2,3,4));A
, , 1

     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6

, , 2

     [,1] [,2] [,3]
[1,]    7    9   11
[2,]    8   10   12

, , 3

     [,1] [,2] [,3]
[1,]   13   15   17
[2,]   14   16   18

, , 4

     [,1] [,2] [,3]
[1,]   19   21   23
[2,]   20   22   24

> B<-aperm(A,c(2,3,1));B
, , 1

     [,1] [,2] [,3] [,4]
[1,]    1    7   13   19
[2,]    3    9   15   21
[3,]    5   11   17   23

, , 2

     [,1] [,2] [,3] [,4]
[1,]    2    8   14   20
[2,]    4   10   16   22
[3,]    6   12   18   24
结果是B把A的第2维移到了第1维,A的第3维移到了第2维,A的第1维移到了第3维,这时有B[i,j,k]=A[i,k,i]
> B[1,3,2];A[3,2,1]
[1] 14
Error in A[3, 2, 1] : 下标出界

6.apply函数

对于向量,可以用sum,mean等函数对其进行计算,对于数组(矩阵),如果想对其1维(或若干维)进行某种计算,可用apply函数,其一般形式为
apply(A,MARGIN,FUN,...)
其中A为一个数组,MARGIN是固定哪些维不变,FUN是用来计算的函数
> A<-matrix(1:6,nrow=2);A
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6
> apply(A, 1, sum)
[1]  9 12
> apply(A, 2, mean)
[1] 1.5 3.5 5.5