一篇关于R矩阵操作总结的神文（终结篇）

2 矩阵计算

2.1矩阵转置

R中矩阵的转置可以用t()函数完成，例如：

> X=matrix(1:12,nrow=4,ncol=3)

> 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

2.2矩阵的行和与列和及行平均值和列均值

在R中很容易计算一个矩阵的各行和和各列和以及各行的平均值和各列的平均值。例如：

> A=matrix(1:12,3,4)

> A

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

[1,]    1    4    7   10

[2,]    2    5    8   11

[3,]    3    6    9   12

> rowSums(A)

[1] 22 26 30

> rowMeans(A)

[1] 5.5 6.5 7.5

> colSums(A)

[1]  6 15 24 33

> colMeans(A)

[1]  2  5  8 11

2.3行列式的值

R中的函数det()将计算方阵A的行列式。例如：

> X=matrix(rnorm(9),nrow=3,ncol=3)

> X

            [,1]       [,2]       [,3]

[1,]  0.05810412 -1.2992698  0.5630315

[2,] -0.28070583  0.1958623 -1.8202283

[3,]  0.83691209  0.4411497  1.0014306

> det(X)

[1] 1.510076

2.4矩阵相加减

矩阵元素的相加减是指维数相同的矩阵，处于同行和同列的位置的元素进行加减。实现这个功能用“＋”，“－”即可。例如：

> A=B=matrix(1:12,nrow=3,ncol=4)

> 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

2.5矩阵的数乘

矩阵的数乘是指一个常数与一个矩阵相乘。设A为m×n矩阵，c≠0，在R中求cA的值，可以用符号“*”。例如：

> c=2

> A=matrix(1:12,nrow=3,ncol=4)

> A

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

[1,]    1    4    7   10

[2,]    2    5    8   11

[3,]    3    6    9   12

> c*A

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

[1,]    2    8   14   20

[2,]    4   10   16   22

[3,]    6   12   18   24

结果矩阵与原矩阵的所有相应元素都差一个常数c。

2.6矩阵相乘

2.6.1矩阵的乘法

A为m×n矩阵，B为n×k矩阵，在R中求AB，可以符号“%*%”。例如：

> A=matrix(1:12,nrow=3,ncol=4)

> B=matrix(1:12,nrow=4,ncol=3)

> A%*%B

     [,1] [,2] [,3]

[1,]   70  158  246

[2,]   80  184  288

[3,]   90  210  330

注意BA无意义，因其不符合矩阵的相乘规则。

若A为n×m矩阵，B为n×k矩阵，在R中求A’B：

> A=matrix(1:12,nrow=4,ncol=3)

> B=matrix(1:12,nrow=4,ncol=3)

> t(A)%*%B

     [,1] [,2] [,3]

[1,]   30   70  110

[2,]   70  174  278

[3,]  110  278  446

也可以用函数crossprod()计算A’B：

> crossprod(A,B)

     [,1] [,2] [,3]

[1,]   30   70  110

[2,]   70  174  278

[3,]  110  278  446


2.6.2矩阵的Kronecker积
 

n×m矩阵A和h×k矩阵B的Kronecker积是一个nh×mk维矩阵，公式为：

              a­11B … a1nB

Am×n×Bh×k=     …     …

              am1B … amnB   mh×nk

在R中Kronecker积可以用函数kronecher()来计算。例如：

> A=matrix(1:4,2,2)

> A

     [,1] [,2]

[1,]    1    3

[2,]    2    4

> B=matrix(rep(1,4),2,2)

> B

     [,1] [,2]

[1,]    1    1

[2,]    1    1

> kronecker(A,B)

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

[1,]    1    1    3    3

[2,]    1    1    3    3

[3,]    2    2    4    4

[4,]    2    2    4    4


2.7矩阵的伴随矩阵

求矩阵A的伴随矩阵可以用LoopAnalyst包中的函数make.adjoint()函数。例如：

>install.packages("LoopAnalyst")

> A=matrix(1:12,nrow=3,ncol=4)

> A

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

[1,]    1    4    7   10

[2,]    2    5    8   11

[3,]    3    6    9   12

> make.adjoint(A)

     [,1] [,2] [,3]

[1,]   -3    6   -3

[2,]    6  -12    6

[3,]   -3    6   -3

2.8矩阵的逆和广义逆

2.8.1矩阵的逆

矩阵A的逆A-1可以用函数solve()，例如：

> A=matrix(rnorm(9),nrow=3,ncol=3)

> A

           [,1]       [,2]        [,3]

[1,] -0.2915845  0.2831544  0.94493154

[2,] -1.6494678  0.6999185 -0.06292334

[3,] -0.7224015 -0.3906971  0.44799963

> solve(A)

          [,1]       [,2]       [,3]

[1,] 0.2359821 -0.4050650 -0.5546321

[2,] 0.6405592  0.4507583 -1.2877720

[3,] 0.9391490 -0.2600663  0.2147417

验证AA-1=1：

> A%*%solve(A)

              [,1]         [,2]          [,3]

[1,]  1.000000e+00 8.433738e-17 -1.341700e-18

[2,]  1.216339e-17 1.000000e+00 -4.667152e-17

[3,] -2.203641e-17 4.283954e-17  1.000000e+00

 

用round函数可以更好的得到结果：

> round(A%*%solve(A))

     [,1] [,2] [,3]

[1,]    1    0    0

[2,]    0    1    0

[3,]    0    0    1

solve()函数也可以用来求解方程组ax=b。

2.8.2矩阵的广义逆（Moore-Penrose）

并非所有的矩阵都有逆，但是所有的矩阵都可有广义逆。n×m矩阵A+是矩阵A的Moore-Penrose逆，如果它满足下列条件：

AA+A=A

A+AA+=A+

(AA+)T=AA+

(A+A)T=A+A

R中MASS包中的ginv()函数可以计算矩阵的Moore-Penrose逆。例如：

> library(MASS)

> A=matrix(1:12,nrow=3,ncol=4)

> A

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

[1,]    1    4    7   10

[2,]    2    5    8   11

[3,]    3    6    9   12

> solve(A)

Error in solve.default(A) : only square matrices can be inverted

> ginv(A)

             [,1]        [,2]        [,3]

[1,] -0.483333333 -0.03333333  0.41666667

[2,] -0.244444444 -0.01111111  0.22222222

[3,] -0.005555556  0.01111111  0.02777778

[4,]  0.233333333  0.03333333 -0.16666667

验证性质①：

> A%*%ginv(A)%*%A

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

[1,]    1    4    7   10

[2,]    2    5    8   11

[3,]    3    6    9   12

> A

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

[1,]    1    4    7   10

[2,]    2    5    8   11

[3,]    3    6    9   12

验证性质②：

>ginv(A)%*%A%*%ginv(A)

           [,1]        [,2]        [,3]

[1,] -0.483333333 -0.03333333  0.41666667

[2,] -0.244444444 -0.01111111  0.22222222

[3,] -0.005555556  0.01111111  0.02777778

[4,]  0.233333333  0.03333333 -0.16666667

> ginv(A)

             [,1]        [,2]        [,3]

[1,] -0.483333333 -0.03333333  0.41666667

[2,] -0.244444444 -0.01111111  0.22222222

[3,] -0.005555556  0.01111111  0.02777778

[4,]  0.233333333  0.03333333 -0.16666667

验证性质③：

> A%*%ginv(A)

          [,1]      [,2]       [,3]

[1,]  0.8333333 0.3333333 -0.1666667

[2,]  0.3333333 0.3333333  0.3333333

[3,] -0.1666667 0.3333333  0.8333333

> t(A%*%ginv(A))

           [,1]      [,2]       [,3]

[1,]  0.8333333 0.3333333 -0.1666667

[2,]  0.3333333 0.3333333  0.3333333

[3,] -0.1666667 0.3333333  0.8333333

验证性质④：

> ginv(A)%*%A

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

[1,]  0.7  0.4  0.1 -0.2

[2,]  0.4  0.3  0.2  0.1

[3,]  0.1  0.2  0.3  0.4

[4,] -0.2  0.1  0.4  0.7

> t(ginv(A)%*%A)

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

[1,]  0.7  0.4  0.1 -0.2

[2,]  0.4  0.3  0.2  0.1

[3,]  0.1  0.2  0.3  0.4

[4,] -0.2  0.1  0.4  0.7

也可以不必如此麻烦来验证性质③和④，因为③和④只是表明AA+和A+A是对称矩阵。

2.8.3 X’X的逆

很多时候，我们需要计算形如X’X的逆。这很容易实现，例如：

> x=matrix(rnorm(9),ncol=3,nrow=3)

> x

           [,1]        [,2]        [,3]

[1,] -0.1806586 -0.76340512 0.002652331

[2,] -1.8018584  0.04467943 1.416332187

[3,]  1.2785359 -1.31653513 0.180653002

> solve(crossprod(x))

          [,1]      [,2]     [,3]

[1,] 1.2181837 0.9664576 1.470940

[2,] 0.9664576 1.2010110 1.204599

[3,] 1.4709402 1.2045986 2.269921

R中的strucchange包中的函数solveCrossprod()也可完成：

> args(solveCrossprod)

function (X, method = c("qr", "chol", "solve"))

NULL

> solveCrossprod(x,method="qr")

         [,1]      [,2]     [,3]

[1,] 1.2181837 0.9664576 1.470940

[2,] 0.9664576 1.2010110 1.204599

[3,] 1.4709402 1.2045986 2.269921

> solveCrossprod(x,method="chol")

          [,1]      [,2]     [,3]

[1,] 1.2181837 0.9664576 1.470940

[2,] 0.9664576 1.2010110 1.204599

[3,] 1.4709402 1.2045986 2.269921

> solveCrossprod(x,method="solve")

          [,1]      [,2]     [,3]

[1,] 1.2181837 0.9664576 1.470940

[2,] 0.9664576 1.2010110 1.204599

[3,] 1.4709402 1.2045986 2.269921

2.9矩阵的特征值和特征向量

可以通过对矩阵A进行谱分解来得到矩阵的特征值和特征向量。矩阵A的谱分解如下：A=UΛU’，其中U的列为A的特征值所对应的特征向量，在R中可以用eigen()函数得到U和Λ。例如：

> args(eigen)

function (x, symmetric, only.values = FALSE, EISPACK = FALSE)

其中，x参数输入矩阵；symmetric参数判断矩阵是否为对称矩阵，如果参数为空，系统将自动检测矩阵的对称性。例如：

> A=matrix(1:9,nrow=3,ncol=3)

> A

     [,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 -4.054214e-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的特征值：

> Aeigen$values

[1]  1.611684e+01 -1.116844e+00 -4.054214e-16

得到矩阵A的特征向量：

> Aeigen$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