R语言里的矩阵处理学习笔记

关于矩阵，通常都会使用matlab来做处理。其实使用R也可以对矩阵做出一些简单的处理。而R语言中提供的matrix，matlab包也提供了不少关于矩阵处理的东西(可以通过??matlab来查看具体函数）。

一、矩阵的输入

通常我们使用函数matrix来创建矩阵，函数的介绍如下：

matrix(data = NA,nrow = 1, ncol = 1, byrow = FALSE, dimnames = NULL)

如果想将数据按行输入矩阵，参数byrow须改为T。

由于矩阵也是特殊的数组，我们也可以用生成数组的函数array（）。具体格式如下：

array(data = NA, dim= length(data), dimnames = NULL)

这里的dim是一个二维数组，生成的就是矩阵了。

当然dim（）函数，attr()（这个是一个格式转换的函数）也是可以生成矩阵的。

还有如diag（）可以产生对角矩阵这样特殊矩阵的函数。

例如生成下面的这个矩阵（为了便于下面的叙述，这个矩阵记为x，生成命令x<-matrix（1:16,2,8））:

[,1]  [,2]  [,3]  [,4]  [,5]  [,6] [,7]  [,8]
[1,]    1    3      5    7    9   11   13   15
[2,]    2    4      6    8  10   12   14   16

我们可以使用如下几种命令：

matrix（1:16,2,8）

x<-1:16  ;dim（x）<-c(2,8)

array(1:16,c(2,8))

x<-1:16;attr(x,"dim")<-c(2,8)

二、矩阵的基本操作

这里主要的有：矩阵的加法与乘法，矩阵求秩，矩阵的转置，方阵求行列式，矩阵求逆，解线性方程组

1、矩阵的加法与乘法

加法直接使用加号，实现对应元素相加。但是要注意两个矩阵必须可加

矩阵的乘法：和matlab类似，R也给出了两种乘法：

“*”：对应位置元素做乘法，如x*x得到结果：

> x*x
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,]    1    9   25   49   81  121 169  225
[2,]    4   16   36   64  100  144  196 256

"%*%":这个是通常意义下的矩阵乘法，如x%*%t(x)得到结果：

> x%*%t(x)
     [,1] [,2]
[1,]  680  744
[2,]  744  816

这里乘法也必须是有意义的才行。

通常我们也使用crossprod()函数来做乘法，crossprod(x,x)效果与x%*%t(x)相同

2、矩阵求秩

这里可以利用qr分解来解决求秩的问题

qr(x)$rank

可以得到x的秩

3、矩阵的转置

常用的命令为t().

R中还可以使用命令aperm()来实现矩阵的广义转置，函数用法如下：

aperm(a, perm, ...)

## Default S3 method:

aperm(a, perm = NULL, resize =TRUE, ...)

## S3 method for class 'table'

aperm(a, perm = NULL, resize =TRUE, keep.class = TRUE, ...)

4、方阵求行列式

命令为det()，无须赘述

5、矩阵求逆与解线性方程组

使用函数solve（）

对于线性方程组b<-A%*%y

的解使用函数solve（A,b)即可

从而我们知道取b为单位阵时即可求逆，通常简化为solve（A)

但是值得注意的是，用直接求逆的办法解线性方程组y<-solve(A)%*%b不仅稳定性低，效率也不咋地。

三、矩阵的分解

这里主要介绍矩阵的lu分解，Cholesky分解，特征值与特征向量，QR分解，奇异值分解

1、LU分解

将系数矩阵A转变成等价两个矩阵L和U的乘积 ，其中L和U分别是下三角和上三角矩阵。当A的所有顺序主子式都不为0时，矩阵A可以分解为A=LU，且分解唯一。其中L是单位下三角矩阵，U是上三角矩阵。

library(Matrix) #这里引入函数包Matrix

> m
     [,1] [,2] [,3]
[1,]    2   -1    3
[2,]    1    2    1
[3,]    2    4    2

> l <- lu(m)

> l
'MatrixFactorization' of Formal class 'denseLU' [package "Matrix"]with 3 slots
  ..@ x   : num [1:9] 2 1 0.5 -1 5 0.5 3 -1 0
  ..@ perm: int [1:3] 1 3 3
  ..@ Dim : int [1:2] 3 3

> LU <- expand(l) #生成P,L,U

> LU
$L
3 x 3 Matrix of class "dtrMatrix" (unitriangular)
     [,1] [,2] [,3]
[1,]  1.0    .    .
[2,]  1.0  1.0    .
[3,]  0.5  0.5  1.0
$U
3 x 3 Matrix of class "dtrMatrix"
     [,1] [,2] [,3]
[1,]    2   -1    3
[2,]    .    5   -1
[3,]    .    .    0
$P（这个矩阵的意思是保持被变换的矩阵第一行不变，二三行对调）
3 x 3 sparse Matrix of class "pMatrix"
          
[1,] | . .
[2,] . . |
[3,] . | .

可以验证A = LU$P%*%LU$L%*%LU$U

P为置换矩阵，L为下单位三角矩阵，U为上三角矩阵；

2、Cholesky分解

如果矩阵A为n阶对称正定矩阵，则存在一个非奇异（满秩）的下三角实矩阵L，使得：A=L%*%t(L)当限定的L的对角元素为正时，分解唯一，成为Cholesky分解

> A
  [,1] [,2] [,3] [,4]
[1,]   2   1   1   1
[2,]   1   2   1   1
[3,]   1   1   2   1
[4,]   1   1   1   2
> chol(A)
        [,1]     [,2]     [,3]    [,4]
[1,] 1.414214 0.7071068 0.7071068 0.7071068
[2,] 0.000000 1.2247449 0.4082483 0.4082483
[3,] 0.000000 0.0000000 1.1547005 0.2886751
[4,] 0.000000 0.0000000 0.0000000 1.1180340
> t(chol(A))%*%chol(A)
  [,1] [,2] [,3] [,4]
[1,]   2   1   1   1
[2,]   1   2   1   1
[3,]   1   1   2   1
[4,]   1   1   1   2
> crossprod(chol(A),chol(A))
  [,1] [,2] [,3] [,4]
[1,]   2   1   1   1
[2,]   1   2   1   1
[3,]   1   1   2   1
[4,]   1   1   1   2
若矩阵为对称正定矩阵，可以利用Cholesky分解求行列式的值，如：
> prod(diag(chol(A))^2)
[1] 5
> det(A)
[1] 5
若矩阵为对称正定矩阵，可以利用Cholesky分解求矩阵的逆，这时用函数chol2inv()，这种用法更有效。

3、特征值与特征向量

函数eigen（A）用来计算方阵A的特征值与特征向量，返回一个含有特征值与特征向量的列表

> A
     [,1] [,2] [,3] [,4]
[1,]    1    5    9   13
[2,]    2    6   10   14
[3,]    3    7   11   15
[4,]    4    8   12   16
> eigen(A)
$values
[1]  3.620937e+01 -2.209373e+00 -1.050249e-15  8.203417e-16
$vectors
           [,1]        [,2]      [,3]       [,4]
[1,] -0.4140028 -0.82289268  0.4422036 -0.1001707
[2,] -0.4688206 -0.42193991 -0.3487083  0.5349238
[3,] -0.5236384 -0.02098714 -0.6291942 -0.7693354
[4,] -0.5784562  0.37996563  0.5356989  0.3345823
有时只需特征值时，我们使用eigen（A,only.value=T)$value可以快速得到结果

4、QR分解

A为m×n矩阵可以进行QR分解，A=QR，其中：Q'Q＝I，在R中可以用函数qr()进行QR分解，例如：
> A=matrix(1:16,4,4) 数据分析培训
> qr(A)
$qr
      [,1]     [,2]       [,3]      [,4]
[1,] -5.4772256 -12.7801930 -2.008316e+01 -2.738613e+01
[2,] 0.3651484 -3.2659863 -6.531973e+00 -9.797959e+00
[3,] 0.5477226 -0.3781696 2.641083e-15 2.056562e-15
[4,] 0.7302967 -0.9124744 8.583032e-01 -2.111449e-16
$rank
[1] 2
$qraux
[1] 1.182574e+00 1.156135e+00 1.513143e+00 2.111449e-16
$pivot
[1] 1 2 3 4
attr(,"class")
[1] "qr"
rank项返回矩阵的秩，qr项包含了矩阵Q和R的信息，要得到矩阵Q和R，可以用函数qr.Q()和qr.R()作用qr()的返回结果，例如：
> qr.R(qr(A))
      [,1]     [,2]       [,3]      [,4]
[1,] -5.477226 -12.780193 -2.008316e+01 -2.738613e+01
[2,] 0.000000 -3.265986 -6.531973e+00 -9.797959e+00
[3,] 0.000000   0.000000 2.641083e-15 2.056562e-15
[4,] 0.000000   0.000000 0.000000e+00 -2.111449e-16
> qr.Q(qr(A))
      [,1]       [,2]     [,3]    [,4]
[1,] -0.1825742 -8.164966e-01 -0.4000874 -0.37407225
[2,] -0.3651484 -4.082483e-01 0.2546329 0.79697056
[3,] -0.5477226 -8.131516e-19 0.6909965 -0.47172438
[4,] -0.7302967 4.082483e-01 -0.5455419 0.04882607
> qr.Q(qr(A))%*%qr.R(qr(A))
  [,1] [,2] [,3] [,4]
[1,]   1   5   9   13
[2,]   2   6   10   14
[3,]   3   7   11   15
[4,]   4   8   12   16

> t(qr.Q(qr(A)))%*%qr.Q(qr(A))
        [,1]       [,2]      [,3]       [,4]
[1,] 1.000000e+00 -1.457168e-16 -6.760001e-17 -7.659550e-17
[2,] -1.457168e-16 1.000000e+00 -4.269046e-17 7.011739e-17
[3,] -6.760001e-17 -4.269046e-17 1.000000e+00 -1.596437e-16
[4,] -7.659550e-17 7.011739e-17 -1.596437e-16 1.000000e+00
> qr.X(qr(A))
  [,1] [,2] [,3] [,4]
[1,]   1   5   9   13
[2,]   2   6   10   14
[3,]   3   7   11   15
[4,]   4   8   12   16

5、svd分解

利用函数svd（）可以对矩阵做svd分解。看一个R提供的例子：

svd> hilbert <- function(n) { i <- 1:n; 1 /outer(i - 1, i, "+") }
svd> X <- hilbert(9)[,1:6]
svd> (s <- svd(X))
$d
[1] 1.668433e+00 2.773727e-01 2.223722e-02 1.084693e-03 3.243788e-05
[6] 5.234864e-07
$u
            [,1]       [,2]       [,3]        [,4]       [,5]        [,6]
 [1,] -0.7244999  0.6265620  0.27350003 -0.08526902 0.02074121 -0.00402455
 [2,] -0.4281556 -0.1298781 -0.64293597  0.55047428 -0.27253421 0.09281592
 [3,] -0.3121985 -0.2803679 -0.33633240 -0.31418014  0.61632113-0.44090375
 [4,] -0.2478932 -0.3141885 -0.06931246 -0.44667149  0.02945426 0.53011986
 [5,] -0.2063780 -0.3140734  0.10786005 -0.30241655 -0.35566839 0.23703838
 [6,] -0.1771408 -0.3026808  0.22105904 -0.09041508 -0.38878613-0.26044927
 [7,] -0.1553452 -0.2877310  0.29280775  0.11551327 -0.19285565-0.42094482
 [8,] -0.1384280 -0.2721599  0.33783778  0.29312535 0.11633231 -0.16079025
 [9,] -0.1248940 -0.2571250  0.36542543  0.43884649 0.46496714  0.43459954
$v
           [,1]       [,2]      [,3]        [,4]       [,5]         [,6]
[1,] -0.7364928  0.6225002  0.2550021 -0.06976287  0.01328234-0.001588146
[2,] -0.4432826 -0.1818705 -0.6866860  0.50860089 -0.19626669  0.041116974
[3,] -0.3274789 -0.3508553 -0.2611139 -0.50473697  0.61605641 -0.259215626
[4,] -0.2626469 -0.3921783  0.1043599 -0.43747940 -0.40833605 0.638901622
[5,] -0.2204199 -0.3945644  0.3509658  0.01612426 -0.46427916-0.675826789
[6,] -0.1904420 -0.3831871  0.5110654  0.53856351  0.44663632 0.257248908
svd> D <- diag(s$d)
svd> s$u %*% D %*% t(s$v) #  X = U D V'
           [,1]      [,2]      [,3]       [,4]       [,5]      [,6]
 [1,] 1.0000000 0.5000000 0.33333333 0.25000000 0.20000000 0.16666667
 [2,] 0.5000000 0.3333333 0.25000000 0.20000000 0.16666667 0.14285714
 [3,] 0.3333333 0.2500000 0.20000000 0.16666667 0.14285714 0.12500000
 [4,] 0.2500000 0.2000000 0.16666667 0.14285714 0.12500000 0.11111111
 [5,] 0.2000000 0.1666667 0.14285714 0.12500000 0.11111111 0.10000000
 [6,] 0.1666667 0.1428571 0.12500000 0.11111111 0.10000000 0.09090909
 [7,] 0.1428571 0.1250000 0.11111111 0.10000000 0.09090909 0.08333333
 [8,] 0.1250000 0.1111111 0.10000000 0.09090909 0.08333333 0.07692308
 [9,] 0.1111111 0.1000000 0.09090909 0.08333333 0.07692308 0.07142857
svd> t(s$u) %*% X %*% s$v #  D = U' X V
              [,1]         [,2]          [,3]         [,4]          [,5]
[1,]  1.668433e+00  2.009230e-16 -2.333610e-16  1.193300e-16 2.347298e-17
[2,]  1.627828e-17  2.773727e-01  7.318365e-19 3.109966e-17 -5.251265e-17
[3,]  1.828617e-17  1.086828e-17  2.223722e-02 4.511721e-18  1.194020e-17
[4,]  2.420517e-17  1.205777e-17  3.433104e-18 1.084693e-03 -4.584063e-18
[5,] -3.406808e-17 -1.147965e-17 -8.968404e-19  6.405788e-18 3.243788e-05
[6,] -1.591696e-17  2.714931e-18  1.721002e-17 -2.358252e-18 1.170640e-17
              [,6]
[1,]  1.015423e-16
[2,]  2.334931e-17
[3,] -1.562373e-17
[4,]  1.364795e-18
[5,] -2.473743e-18
[6,]  5.234864e-07

6、矩阵广义逆(Moore-Penrose)
  n×m矩阵A+称为m×n矩阵A的Moore-Penrose逆，如果它满足下列条件：
①   A A+A=A；②A+A A+= A+；③(A A+)H=A A+；④(A+A)H= A+A
在R的MASS包中的函数ginv()可计算矩阵A的Moore-Penrose逆，例如：
library(“MASS”)
> A
  [,1] [,2] [,3] [,4]
[1,]   1   5   9   13
[2,]   2   6   10   14
[3,]   3   7   11   15
[4,]   4   8   12   16
> ginv(A)
    [,1]   [,2] [,3]   [,4]
[1,] -0.285 -0.1075 0.07 0.2475
[2,] -0.145 -0.0525 0.04 0.1325
[3,] -0.005 0.0025 0.01 0.0175
[4,] 0.135 0.0575 -0.02 -0.0975
验证性质1：
> A%*%ginv(A)%*%A
  [,1] [,2] [,3] [,4]
[1,]   1   5   9   13
[2,]   2   6   10   14
[3,]   3   7   11   15
[4,]   4   8   12   16
验证性质2：
> ginv(A)%*%A%*%ginv(A)
    [,1]   [,2] [,3]   [,4]
[1,] -0.285 -0.1075 0.07 0.2475
[2,] -0.145 -0.0525 0.04 0.1325
[3,] -0.005 0.0025 0.01 0.0175
[4,] 0.135 0.0575 -0.02 -0.0975
验证性质3:
> t(A%*%ginv(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
> A%*%ginv(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
验证性质4:
> 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
> 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

对于矩阵，我们还可以使用cbind(),rbind()构造分块矩阵。