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

3 矩阵高级操作

3.1 Choleskey分解

对于正定矩阵A，可以对其进行Choleskey分解，A=P’P，其中P为上三角矩阵，在R中可以用函数chol()。例如：

> A=diag(3)+1

> A

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

[1,]    2    1    1

[2,]    1    2    1

[3,]    1    1    2

> chol(A)

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

[1,] 1.414214 0.7071068 0.7071068

[2,] 0.000000 1.2247449 0.4082483

[3,] 0.000000 0.0000000 1.1547005

验证A=P’P：

> t(chol(A))%*%chol(A)

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

[1,]    2    1    1

[2,]    1    2    1

[3,]    1    1    2

也可以用crossprod()函数：

> crossprod(chol(A),chol(A))

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

[1,]    2    1    1

[2,]    1    2    1

[3,]    1    1    2

可以用Choleskey分解来计算矩阵的行列式：

> prod(diag(chol(A))^2)

[1] 4

> det(A)

[1] 4

也可以用Choleskey分解来计算矩阵的逆，这时候可以用到函数chol2inv():

> chol2inv(chol(A))

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

[1,]  0.75 -0.25 -0.25

[2,] -0.25  0.75 -0.25

[3,] -0.25 -0.25  0.75

> solve(A)

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

[1,]  0.75 -0.25 -0.25

[2,] -0.25  0.75 -0.25

[3,] -0.25 -0.25  0.75

3.2奇异值分解

A为m×n矩阵，矩阵的秩为r。A可以分解为A=UDV’，其中U’U=V’V=I。在R中可以用函数svd()。例如：

> A=matrix(1:18,3,6)

> A

     [,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)

$d

[1] 4.589453e+01 1.640705e+00 2.294505e-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.7196003 -0.67039144

[2,] -0.19033085 -0.5089325  0.55766549

[3,] -0.30329950 -0.2982646  0.28189237

[4,] -0.41626816 -0.0875968  0.07320847

[5,] -0.52923682  0.1230711  0.12920119

[6,] -0.64220548  0.3337389 -0.37157608

> A.u%*%diag(A.d)%*%t(A.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

3.3 QR分解

A为m×n矩阵可以进行QR分解:A=QR，其中Q’Q=I，在R中可以用函数qr()来完成这个过程，例如：

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

> qr(A)

$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  7.880925e-16

[4,]  0.7302967  -0.9124744  9.277920e-01

$rank

[1] 2

$qraux

[1] 1.182574 1.156135 1.373098

$pivot

[1] 1 2 3

attr(,"class")

[1] "qr"

Rank返回的是矩阵的秩。Qr项包含了Q矩阵和R矩阵的信息。要想得到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  4.938541e-17  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  7.880925e-16