专栏首页优雅RR 语言中的矩阵计算

R 语言中的矩阵计算

这方面有需求

很棒的内容,转载不了我就直接扣过来了~~对这方面有需求的读者不妨收藏一下。

作者:张丹(Conan) 来源:http://blog.fens.me/r-matrix/

前言

R 是作为统计语言,生来就对数学有良好的支持。矩阵计算作为底层的数学工具,有非常广泛的使用场景。用R语言很好地封装了,矩阵的各种计算方法,一个函数一行代码,就能完成复杂的矩阵分解等操作。让建模人员可以更专注于模型推理和业务逻辑实现,把复杂的矩阵计算交给R语言来完成。 本文总结了 R 语言用于矩阵的各种计算操作。

1. 基本操作

# 生成矩阵 
> m<-matrix(1:20,4,5);m
     [,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

# 矩阵行
> nrow(m)
[1] 4
# 矩阵列
> ncol(m)
[1] 5

取对角线元素,生成对角矩阵:

# 对角线元素
> diag(m)
[1]  1  6 11 16

# 以对角线元素,生成对角矩阵
> diag(diag(m))
     [,1] [,2] [,3] [,4]
[1,]    1    0    0    0
[2,]    0    6    0    0
[3,]    0    0   11    0

上三角,下三角:

# 上三角
> m<-matrix(1:20,4,5)
> upper.tri(m)
      [,1]  [,2]  [,3]  [,4] [,5]
[1,] FALSE  TRUE  TRUE  TRUE TRUE
[2,] FALSE FALSE  TRUE  TRUE TRUE
[3,] FALSE FALSE FALSE  TRUE TRUE
[4,] FALSE FALSE FALSE FALSE TRUE

# 上三角值
> m[which(upper.tri(m))] 
 [1]  5  9 10 13 14 15 17 18 19 20

# 上三角矩阵
> m[!upper.tri(m)]<-0;m
     [,1] [,2] [,3] [,4] [,5]
[1,]    0    5    9   13   17
[2,]    0    0   10   14   18
[3,]    0    0    0   15   19
[4,]    0    0    0    0   20

# 下三角矩阵
> m[!lower.tri(m)]<-0;m
     [,1] [,2] [,3] [,4] [,5]
[1,]    0    0    0    0    0
[2,]    2    0    0    0    0
[3,]    3    7    0    0    0
[4,]    4    8   12    0    0

矩阵转置:

> m<-matrix(1:20,4,5);m
     [,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

# 转置,行列互转
> t(m)
     [,1] [,2] [,3] [,4]
[1,]    1    2    3    4
[2,]    5    6    7    8
[3,]    9   10   11   12
[4,]   13   14   15   16
[5,]   17   18   19   20

对角矩阵填充:

# 创建方阵
> m<-matrix(1:16,4,4);m
     [,1] [,2] [,3] [,4]
[1,]    1    5    9   13
[2,]    2    6   10   14
[3,]    3    7   11   15
[4,]    4    8   12   16

# 用上三角填充下三角
> m[lower.tri(m)]<-m[upper.tri(m)]
> m
     [,1] [,2] [,3] [,4]
[1,]    1    5    9   13
[2,]    5    6   10   14
[3,]    9   13   11   15
[4,]   10   14   15   16

填充后,发现矩阵并不是对称的,原因是上三角取值按列取值,所以先取 10 后取 13,导致上三角和下三角取值顺序不完全一致。

> m<-matrix(1:16,4,4);m
     [,1] [,2] [,3] [,4]
[1,]    1    5    9   13
[2,]    2    6   10   14
[3,]    3    7   11   15
[4,]    4    8   12   16

# 上三角取值
> m[upper.tri(m)]
[1]  5  9 10 13 14 15

# 下三角取值
> m[lower.tri(m)]
[1]  2  3  4  7  8 12

调整后,我们要先转置,再取值再填充,形成对称结构。

> m<-matrix(1:20,4,5)

# 转置后,取下三角,填充下三角
> m[lower.tri(m)]<-t(m)[lower.tri(m)]
> m
     [,1] [,2] [,3] [,4]
[1,]    1    5    9   13
[2,]    5    6   10   14
[3,]    9   10   11   15
[4,]   13   14   15   16

矩阵和 data.frame 转换,用行列形成索引结构。

> m<-matrix(1:12,4,3);m
     [,1] [,2] [,3]
[1,]    1    5    9
[2,]    2    6   10
[3,]    3    7   11
[4,]    4    8   12

# 矩阵转行列data.frame
> row<-rep(1:nrow(m),ncol(m))              # 行
> col<-rep(1:ncol(m),each=nrow(m))         # 列
> df<-data.frame(row,col,v=as.numeric(m))  # 行列索引数据框
> df
   row col  v
1    1   1  1
2    2   1  2
3    3   1  3
4    4   1  4
5    1   2  5
6    2   2  6
7    3   2  7
8    4   2  8
9    1   3  9
10   2   3 10
11   3   3 11
12   4   3 12

# 行列索引数据框转矩阵
> m<-matrix(0,length(unique(df$row)),length(unique(df$col)) )
> apply(df,1,function(dat){
+     m[dat[1],dat[2]]<<-dat[3]  
+     invisible()
+ })
# 打印矩阵
> m
     [,1] [,2] [,3]
[1,]    1    5    9
[2,]    2    6   10
[3,]    3    7   11
[4,]    4    8   12

2. 矩阵计算

加法,减法。

# 加载矩阵计算工具包
> library(matrixcalc)

# 新建2个矩阵,行列长度相同
> m0<-matrix(1:20,4,5);m0
     [,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
> m1<-matrix(sample(100,20),4,5);m1
     [,1] [,2] [,3] [,4] [,5]
[1,]   40   79   97   57   78
[2,]   93   32   48    8   95
[3,]   63    6   56   12    9
[4,]   28   31   72   27   26

# 矩阵加法
> m0+m1
     [,1] [,2] [,3] [,4] [,5]
[1,]   41   84  106   70   95
[2,]   95   38   58   22  113
[3,]   66   13   67   27   28
[4,]   32   39   84   43   46

# 矩阵减法
> m0-m1
     [,1] [,2] [,3] [,4] [,5]
[1,]  -39  -74  -88  -44  -61
[2,]  -91  -26  -38    6  -77
[3,]  -60    1  -45    3   10
[4,]  -24  -23  -60  -11   -6

矩阵值相乘。

> m0*m1
     [,1] [,2] [,3] [,4] [,5]
[1,]   40  395  873  741 1326
[2,]  186  192  480  112 1710
[3,]  189   42  616  180  171
[4,]  112  248  864  432  520

矩阵乘法,满足第二个矩阵的列数和第一个矩阵的行数相等,所以把上面生成的 m0 矩阵( 4 行 5 列)转置为( 5 行 4 列),再用 m1 矩阵( 4 行 5 列),进行矩阵乘法,得到一个 5 行 5 列的结果矩阵。

> t(m0)%*%m1
     [,1] [,2] [,3] [,4] [,5]
[1,]  527  285  649  217  399
[2,] 1423  877 1741  633 1231
[3,] 2319 1469 2833 1049 2063
[4,] 3215 2061 3925 1465 2895
[5,] 4111 2653 5017 1881 3727

# 通过函数实现矩阵相乘
> crossprod(m0,m1)
     [,1] [,2] [,3] [,4] [,5]
[1,]  527  285  649  217  399
[2,] 1423  877 1741  633 1231
[3,] 2319 1469 2833 1049 2063
[4,] 3215 2061 3925 1465 2895
[5,] 4111 2653 5017 1881 3727

矩阵外积。

> m<-matrix(1:6,2);m
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6
> t(m)
     [,1] [,2]
[1,]    1    2
[2,]    3    4
[3,]    5    6

> m %o% t(m)  # 外积,同outer(m,t(m))
, , 1, 1

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

, , 2, 1

     [,1] [,2] [,3]
[1,]    3    9   15
[2,]    6   12   18

, , 3, 1

     [,1] [,2] [,3]
[1,]    5   15   25
[2,]   10   20   30

, , 1, 2

     [,1] [,2] [,3]
[1,]    2    6   10
[2,]    4    8   12

, , 2, 2

     [,1] [,2] [,3]
[1,]    4   12   20
[2,]    8   16   24

, , 3, 2

     [,1] [,2] [,3]
[1,]    6   18   30
[2,]   12   24   36

矩阵直和。

> m0<-matrix(1:4,2,2);m0
     [,1] [,2]
[1,]    1    3
[2,]    2    4
> m1<-matrix(1:9,3,3);m1
     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9

> m0 %s% m1       #矩阵直和,同 direct.sum(m0,m1) 
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    3    0    0    0
[2,]    2    4    0    0    0
[3,]    0    0    1    4    7
[4,]    0    0    2    5    8
[5,]    0    0    3    6    9

矩阵直积。

> m0<-matrix(1:4,2,2);m0
     [,1] [,2]
[1,]    1    3
[2,]    2    4
> m1<-matrix(1:9,3,3);m1
     [,1] [,2] [,3]
[1,]    1    4    7
[2,]    2    5    8
[3,]    3    6    9

> m0 %x% m1     #矩阵直积,同 kronecker(m0,m1),direct.prod(m0,m1)
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    4    7    3   12   21
[2,]    2    5    8    6   15   24
[3,]    3    6    9    9   18   27
[4,]    2    8   14    4   16   28
[5,]    4   10   16    8   20   32
[6,]    6   12   18   12   24   36

3. 矩阵性质

3.1 奇异矩阵

首先,我们线创建一个非奇异矩阵,判断非奇异矩阵方法,行列式不等于 0,矩阵可逆,满秩。

# 创建一个非奇异矩阵
> m1<-matrix(sample(1:16),4)

# 行列式不等于0
> det(m1)     # !=0
[1] 14193

# 有逆矩阵
> solve(m1)   # 可逆
             [,1]        [,2]        [,3]         [,4]
[1,] -0.026210104  0.09666737  0.02473050 -0.119988727
[2,] -0.002325090  0.08384415 -0.07038681  0.008173043
[3,] -0.007186641 -0.13478475  0.05516804  0.176777285
[4,]  0.074614246  0.03663778 -0.01395054 -0.080462200

# 满秩
> library(matrixcalc)
> matrix.rank(m1)    # 长度=n,满秩
[1] 4

再创建一个奇异矩阵,判断奇异矩阵方法包括,行列式等于 0,矩阵不可逆,不是满秩。

# 奇异矩阵
> m0<-matrix(1:16,4)

# 计算行列式
> det(m0)     
[1] 0

# 不可逆
> solve(m0)   
Error in solve.default(m0) : 
  Lapack routine dgesv: system is exactly singular: U[3,3] = 0

# 不是满秩
> matrix.rank(m0)     # 长度<4
[1] 2

3.2 逆矩阵

# 创建方阵,非奇异矩阵
> m0<-matrix(sample(100,16),4,4);m0
     [,1] [,2] [,3] [,4]
[1,]   24   31   80   37
[2,]   84   13   42   71
[3,]   95   62   93   86
[4,]   69   16   94   35

# 计算矩阵的逆矩阵
> solve(m0)
            [,1]          [,2]        [,3]         [,4]
[1,] -0.03083680 -0.0076561475  0.01258023  0.017218514
[2,] -0.01710957 -0.0270246488  0.03152548 -0.004553923
[3,]  0.01384721 -0.0003070371 -0.00886117  0.007757524
[4,]  0.03142440  0.0282722871 -0.01541411 -0.024126340

# 逆矩阵的性质,逆矩阵与原矩阵进行矩阵乘法,得到对角矩阵。
> round(solve(m0) %*% m0)  # 对角矩阵
     [,1] [,2] [,3] [,4]
[1,]    1    0    0    0
[2,]    0    1    0    0
[3,]    0    0    1    0
[4,]    0    0    0    1

广义逆矩阵,将逆矩阵的概率推广到奇异矩阵和长方形矩阵上,就产生了广义逆矩阵。

# 创建奇异矩阵
> m<-matrix(1:16,4,4);m 
     [,1] [,2] [,3] [,4]
[1,]    1    5    9   13
[2,]    2    6   10   14
[3,]    3    7   11   15
[4,]    4    8   12   16

# 逆矩阵计算失败
> solve(m)
Error in solve.default(m) : 
  Lapack routine dgesv: system is exactly singular: U[3,3] = 0

# 广义逆矩阵
> library(MASS) # 加载MASS包
> ginv(m)
       [,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

ginv 函数计算非奇异矩阵,和 solve() 函数比较。

# 非奇异矩阵
> m0<-matrix(sample(100,16),4,4);m0
     [,1] [,2] [,3] [,4]
[1,]    5   61   75   74
[2,]   10   67   11   21
[3,]   29   32   92   17
[4,]   72   23   87   36

# 计算广义矩阵的逆矩阵,与solve()结果相同
> ginv(m0)
              [,1]         [,2]         [,3]         [,4]
[1,] -0.0080590817  0.006534814 -0.011333076  0.018105645
[2,] -0.0041284695  0.017615769  0.005333304 -0.004308072
[3,]  0.0009226026 -0.006671759  0.016312821 -0.005707878
[4,]  0.0165261737 -0.008200731 -0.020163889  0.008112906

# 计算矩阵的逆矩阵
> solve(m0)
              [,1]         [,2]         [,3]         [,4]
[1,] -0.0080590817  0.006534814 -0.011333076  0.018105645
[2,] -0.0041284695  0.017615769  0.005333304 -0.004308072
[3,]  0.0009226026 -0.006671759  0.016312821 -0.005707878
[4,]  0.0165261737 -0.008200731 -0.020163889  0.008112906

逆矩阵的特性。

> m<-matrix(1:16,4,4);m 
     [,1] [,2] [,3] [,4]
[1,]    1    5    9   13
[2,]    2    6   10   14
[3,]    3    7   11   15
[4,]    4    8   12   16

# 原矩阵%*%广义逆矩阵%*%原矩阵=原矩阵
> m %*% ginv(m) %*% m 
     [,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(m) %*% m %*% ginv(m) 
       [,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

# 原矩阵%*%广义逆矩阵=转置后的,原矩阵%*%广义逆矩阵
> m %*% ginv(m)    
     [,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(m %*% ginv(m)) 
     [,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(m) %*% m   
     [,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(m) %*% m)
     [,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

3.3 特征值和特征向量

# 创建一个方阵
> m0<-matrix(sample(100,16),4,4);m0
     [,1] [,2] [,3] [,4]
[1,]   97   93   11   70
[2,]   19   58   23   90
[3,]   17   94    6   35
[4,]   79   71   43   88

# 计算特征值和特征向量
> eigen(m0)
eigen() decomposition
$values
[1] 236.14449+ 0.00000i  40.51501+ 0.00000i -13.82975+32.81166i -13.82975-32.81166i

$vectors
             [,1]          [,2]                  [,3]                  [,4]
[1,] 0.6055193+0i -0.6428709+0i  0.2114869-0.0613776i  0.2114869+0.0613776i
[2,] 0.4115920+0i  0.3939259+0i -0.0781518+0.3993580i -0.0781518-0.3993580i
[3,] 0.3054253+0i  0.6482306+0i  0.7557355+0.0000000i  0.7557355+0.0000000i
[4,] 0.6088134+0i -0.1064728+0i -0.3210016-0.3342655i -0.3210016+0.3342655i

symmetric=TRUE 时,计算对称矩阵的特征值和特征向量,当 m0 不是对称矩阵时,则取下三角对称结构进行计算。

> eigen(m0,symmetric = TRUE)
eigen() decomposition
$values
[1] 231.824838  87.176816  -3.422899 -66.578756

$vectors
           [,1]       [,2]       [,3]        [,4]
[1,] -0.4801643  0.7155875  0.5040624 -0.05742733
[2,] -0.5024315 -0.5470752  0.2262475 -0.63014551
[3,] -0.3634220 -0.4138943  0.3287896  0.76714627
[4,] -0.6204267  0.1316618 -0.7659181  0.10538188

# 生成下三角矩阵的对称矩阵
> m0[upper.tri(m0)]<-t(m0)[upper.tri(m0)];m0
     [,1] [,2] [,3] [,4]
[1,]   97   19   17   79
[2,]   19   58   94   71
[3,]   17   94    6   43
[4,]   79   71   43   88

# 计算特征值,与eigen(m0,symmetric = TRUE) 一致
> eigen(m0)
eigen() decomposition
$values
[1] 231.824838  87.176816  -3.422899 -66.578756

$vectors
           [,1]       [,2]       [,3]        [,4]
[1,] -0.4801643  0.7155875  0.5040624 -0.05742733
[2,] -0.5024315 -0.5470752  0.2262475 -0.63014551
[3,] -0.3634220 -0.4138943  0.3287896  0.76714627
[4,] -0.6204267  0.1316618 -0.7659181  0.10538188

4. 矩阵分解

下面将介绍 4 种矩阵常用的分解的方法,包括三角分解 LU,choleskey 分解,QR 分解,奇异值分解 SVD。

4.1 三角分解 LU

三角分解法是将原方阵分解成一个上三角形矩阵和一个下三角形矩阵,这样的分解法又称为 LU 分解法。它的用途主要在简化一个大矩阵的行列式值的计算过程,求逆矩阵,和求解联立方程组。这种分解法所得到的上下三角形矩阵不唯一,一对上三角形矩阵和下三角形矩阵,矩阵相乘会得到原矩阵。

# 创建一个矩阵
> m0<-matrix(sample(100,16),4);m0
     [,1] [,2] [,3] [,4]
[1,]   25   88   35   87
[2,]   26   75   73   15
[3,]   36   62   80   53
[4,]   70   44   21   47

# 三角分解
> lu<-lu.decomposition(m0)
> lu
$L                                        # 下三角矩阵
     [,1]      [,2]     [,3] [,4]
[1,] 1.00  0.000000 0.000000    0
[2,] 1.04  1.000000 0.000000    0
[3,] 1.44  3.917676 1.000000    0
[4,] 2.80 12.251816 4.617547    1

$U                                        # 上三角矩阵
     [,1]   [,2]      [,3]      [,4]
[1,]   25  88.00   35.0000   87.0000
[2,]    0 -16.52   36.6000  -75.4800
[3,]    0   0.00 -113.7869  223.4262
[4,]    0   0.00    0.0000 -303.5137

# 三角分解的下三角矩阵L %*% 三角分解的上三角矩阵U = 原矩阵
> lu$L %*% lu$U
     [,1] [,2] [,3] [,4]
[1,]   25   88   35   87
[2,]   26   75   73   15
[3,]   36   62   80   53
[4,]   70   44   21   47

4.2 choleskey 分解

Cholesky 分解是把一个对称正定的矩阵表示成一个下三角矩阵L和其转置的乘积的分解。它要求矩阵的所有特征值必须大于零,故分解的下三角的对角元也是大于零的。Cholesky 分解法又称平方根法,是当A为实对称正定矩阵时,LU 三角分解法的变形。

# 创建对称方阵
> m1<-diag(5)+2;m1
     [,1] [,2] [,3] [,4] [,5]
[1,]    3    2    2    2    2
[2,]    2    3    2    2    2
[3,]    2    2    3    2    2
[4,]    2    2    2    3    2
[5,]    2    2    2    2    3

# choleskey分解
> chol(m1)
         [,1]     [,2]      [,3]      [,4]      [,5]
[1,] 1.732051 1.154701 1.1547005 1.1547005 1.1547005
[2,] 0.000000 1.290994 0.5163978 0.5163978 0.5163978
[3,] 0.000000 0.000000 1.1832160 0.3380617 0.3380617
[4,] 0.000000 0.000000 0.0000000 1.1338934 0.2519763
[5,] 0.000000 0.000000 0.0000000 0.0000000 1.1055416

# choleskey分解的性质
# chol分解的逆矩阵 %*% chol分解矩阵 = 原矩阵
> t(chol(m1))%*%chol(m1)
     [,1] [,2] [,3] [,4] [,5]
[1,]    3    2    2    2    2
[2,]    2    3    2    2    2
[3,]    2    2    3    2    2
[4,]    2    2    2    3    2
[5,]    2    2    2    2    3

# chol分解矩阵的对角线值的平方的积 = 行列式的模数
> prod(diag(chol(m1))^2)
[1] 11
> det(m1)
[1] 11

# chol分解的逆
> chol2inv(m1)
             [,1]         [,2]        [,3]         [,4]         [,5]
[1,]  0.166658199 -0.055580958 -0.01859473 -0.006401463 -0.002743484
[2,] -0.055580958  0.166590459 -0.05578418 -0.019204390 -0.008230453
[3,] -0.018594726 -0.055784179  0.16598080 -0.057613169 -0.024691358
[4,] -0.006401463 -0.019204390 -0.05761317  0.160493827 -0.074074074
[5,] -0.002743484 -0.008230453 -0.02469136 -0.074074074  0.111111111

> chol2inv(chol(m1))
           [,1]       [,2]       [,3]       [,4]       [,5]
[1,]  0.8181818 -0.1818182 -0.1818182 -0.1818182 -0.1818182
[2,] -0.1818182  0.8181818 -0.1818182 -0.1818182 -0.1818182
[3,] -0.1818182 -0.1818182  0.8181818 -0.1818182 -0.1818182
[4,] -0.1818182 -0.1818182 -0.1818182  0.8181818 -0.1818182
[5,] -0.1818182 -0.1818182 -0.1818182 -0.1818182  0.8181818

4.3 QR 分解

QR 分解法是将矩阵分解成一个正规正交矩阵与上三角形矩阵,所以称为 QR 分解法,与此正规正交矩阵的通用符号 Q 有关。

# 创建对称方阵
> m1<-diag(5)+2;m1
     [,1] [,2] [,3] [,4] [,5]
[1,]    3    2    2    2    2
[2,]    2    3    2    2    2
[3,]    2    2    3    2    2
[4,]    2    2    2    3    2
[5,]    2    2    2    2    3

# QR分解
> q<-qr(m1);q
$qr
     [,1]       [,2]       [,3]       [,4]       [,5]
[1,] -5.0 -4.8000000 -4.8000000 -4.8000000 -4.8000000
[2,]  0.4 -1.4000000 -0.6857143 -0.6857143 -0.6857143
[3,]  0.4  0.2142857 -1.2205720 -0.4012839 -0.4012839
[4,]  0.4  0.2142857  0.1560549 -1.1527216 -0.2852095
[5,]  0.4  0.2142857  0.1560549  0.1246843  1.1168808

$rank
[1] 5

$qraux
[1] 1.600000 1.928571 1.975343 1.992196 1.116881

$pivot
[1] 1 2 3 4 5

attr(,"class")
[1] "qr"

# 计算QR分解矩阵的R值
> qr.R(q)
     [,1] [,2]       [,3]       [,4]       [,5]
[1,]   -5 -4.8 -4.8000000 -4.8000000 -4.8000000
[2,]    0 -1.4 -0.6857143 -0.6857143 -0.6857143
[3,]    0  0.0 -1.2205720 -0.4012839 -0.4012839
[4,]    0  0.0  0.0000000 -1.1527216 -0.2852095
[5,]    0  0.0  0.0000000  0.0000000  1.1168808

# 计算QR分解矩阵的Q值
> qr.Q(q)
     [,1]        [,2]        [,3]        [,4]       [,5]
[1,] -0.6  0.62857143  0.36784361  0.26144202 -0.2030692
[2,] -0.4 -0.77142857  0.36784361  0.26144202 -0.2030692
[3,] -0.4 -0.05714286 -0.85272836  0.26144202 -0.2030692
[4,] -0.4 -0.05714286 -0.03344033 -0.89127960 -0.2030692
[5,] -0.4 -0.05714286 -0.03344033 -0.02376746  0.9138115

# QR分解的性质
# Q的值 %*% R值 = 原矩阵
> qr.Q(q) %*% qr.R(q) #=m1
     [,1] [,2] [,3] [,4] [,5]
[1,]    3    2    2    2    2
[2,]    2    3    2    2    2
[3,]    2    2    3    2    2
[4,]    2    2    2    3    2
[5,]    2    2    2    2    3

# X值 = 原矩阵,同Q的值 %*% R值 
> qr.X(q) #=m1
     [,1] [,2] [,3] [,4] [,5]
[1,]    3    2    2    2    2
[2,]    2    3    2    2    2
[3,]    2    2    3    2    2
[4,]    2    2    2    3    2
[5,]    2    2    2    2    3

# Q值的转置矩阵 * Q值 = 
> t(qr.Q(q)) %*% qr.Q(q)
              [,1]          [,2]          [,3]         [,4]          [,5]
[1,]  1.000000e+00 -3.469447e-17 -2.081668e-17 1.214306e-17  5.551115e-17
[2,] -3.469447e-17  1.000000e+00  7.199102e-17 5.204170e-17 -7.632783e-17
[3,] -2.081668e-17  7.199102e-17  1.000000e+00 3.382711e-17 -6.938894e-18
[4,]  1.214306e-17  5.204170e-17  3.382711e-17 1.000000e+00  3.122502e-17
[5,]  5.551115e-17 -7.632783e-17 -6.938894e-18 3.122502e-17  1.000000e+00

4.4 奇异值分解 SVD

奇异值分解 (singular value decomposition, SVD) 是一种正交矩阵分解法。SVD 是最可靠的分解法,但是它比 QR 分解法要花上近十倍的计算时间。[U,S,V]=svd(A),其中 UV 分别代表两个正交矩阵,而 S 代表一对角矩阵。和 QR 分解法相同, 原矩阵 A 不必为正方矩阵。使用 SVD 分解法的用途是解最小平方误差法和数据压缩。

# 创建对称方阵
> m1<-diag(5)+2;m1
     [,1] [,2] [,3] [,4] [,5]
[1,]    3    2    2    2    2
[2,]    2    3    2    2    2
[3,]    2    2    3    2    2
[4,]    2    2    2    3    2
[5,]    2    2    2    2    3

# 奇异值分解
> s<-svd(m1);s
$d
[1] 11  1  1  1  1

$u
           [,1]          [,2]          [,3]          [,4]       [,5]
[1,] -0.4472136  4.440892e-17 -3.330669e-17 -3.108624e-16  0.8944272
[2,] -0.4472136 -3.219647e-16  2.414735e-16  8.660254e-01 -0.2236068
[3,] -0.4472136 -5.773503e-01  5.773503e-01 -2.886751e-01 -0.2236068
[4,] -0.4472136 -2.113249e-01 -7.886751e-01 -2.886751e-01 -0.2236068
[5,] -0.4472136  7.886751e-01  2.113249e-01 -2.886751e-01 -0.2236068

$v
           [,1]          [,2]          [,3]       [,4]       [,5]
[1,] -0.4472136  0.000000e+00  0.000000e+00  0.0000000  0.8944272
[2,] -0.4472136 -1.665335e-16  1.457168e-16  0.8660254 -0.2236068
[3,] -0.4472136 -5.773503e-01  5.773503e-01 -0.2886751 -0.2236068
[4,] -0.4472136 -2.113249e-01 -7.886751e-01 -0.2886751 -0.2236068
[5,] -0.4472136  7.886751e-01  2.113249e-01 -0.2886751 -0.2236068

# 奇异分解性质
# svd的u矩阵 %*% svd的d矩阵的对角性值 %*% svd的v的转置矩阵 = 原矩阵
> s$u %*% diag(s$d) %*% t(s$v) #=m1
[,1] [,2] [,3] [,4] [,5]
[1,]    3    2    2    2    2
[2,]    2    3    2    2    2
[3,]    2    2    3    2    2
[4,]    2    2    2    3    2
[5,]    2    2    2    2    3

# svd的u矩阵的转置矩阵 %*% svd的u矩阵 = svd的s矩阵的转置矩阵 %*% svd的v矩阵
> t(s$u) %*% s$u
[,1]          [,2]          [,3]          [,4]          [,5]
[1,]  1.000000e+00  0.000000e+00 -5.551115e-17  0.000000e+00 -2.775558e-17
[2,]  0.000000e+00  1.000000e+00 -2.775558e-17 -1.387779e-16  8.326673e-17
[3,] -5.551115e-17 -2.775558e-17  1.000000e+00  6.245005e-17 -6.938894e-17
[4,]  0.000000e+00 -1.387779e-16  6.245005e-17  1.000000e+00 -1.387779e-16
[5,] -2.775558e-17  8.326673e-17 -6.938894e-17 -1.387779e-16  1.000000e+00
> t(s$v) %*% s$v
[,1]          [,2]         [,3]          [,4]          [,5]
[1,]  1.000000e+00  0.000000e+00 5.551115e-17 -1.665335e-16  2.775558e-17
[2,]  0.000000e+00  1.000000e+00 8.326673e-17 -8.326673e-17  0.000000e+00
[3,]  5.551115e-17  8.326673e-17 1.000000e+00  1.595946e-16  2.775558e-17
[4,] -1.665335e-16 -8.326673e-17 1.595946e-16  1.000000e+00 -8.326673e-17
[5,]  2.775558e-17  0.000000e+00 2.775558e-17 -8.326673e-17  1.000000e+00

5. 特殊矩阵

下面介绍的多种特殊矩阵,都是在 matrixcalc 库中提供的。

5.1 Hankel Matrix

汉克尔矩阵 (Hankel Matrix) 是具有恒定倾斜对角线的方形矩阵。Hankel 矩阵的行列式称为 catalecticant。该函数根据 n 向量 x 的值构造 n 阶 Hankel 矩阵。矩阵的每一行是前一行中值的循环移位。 矩阵定义:

> hankel.matrix(6, 1:50)
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    2    3    4    5    6
[2,]    2    3    4    5    6    7
[3,]    3    4    5    6    7    8
[4,]    4    5    6    7    8    9
[5,]    5    6    7    8    9   10
[6,]    6    7    8    9   10   11

5.2 Hilbert Matrix

希尔伯特矩阵是一种数学变换矩阵,正定,且高度病态(即,任何一个元素发生一点变动,整个矩阵的行列式的值和逆矩阵都会发生巨大变化),病态程度和阶数相关。希尔伯特矩阵是一种特殊的汉克尔矩阵,该函数返回 n 乘 n 希尔伯特矩阵。 矩阵定义:

> hilbert.matrix(4)
          [,1]      [,2]      [,3]      [,4]
[1,] 1.0000000 0.5000000 0.3333333 0.2500000
[2,] 0.5000000 0.3333333 0.2500000 0.2000000
[3,] 0.3333333 0.2500000 0.2000000 0.1666667
[4,] 0.2500000 0.2000000 0.1666667 0.1428571

5.3 Creation Matrix

创造矩阵,n 阶创建矩阵也称为推导矩阵,该函数返回阶数 n 创建矩阵,在主对角线下方的子对角线上具有序列 1,2,…,n-1 的方阵。 矩阵定义:

> creation.matrix(5)  
     [,1] [,2] [,3] [,4] [,5]
[1,]    0    0    0    0    0
[2,]    1    0    0    0    0
[3,]    0    2    0    0    0
[4,]    0    0    3    0    0
[5,]    0    0    0    4    0

5.4 Stirling Matrix

斯特林公式(Stirling’s approximation)是一条用来取n的阶乘的近似值的数学公式。一般来说,当 n 很大的时候,n 阶乘的计算量十分大,所以斯特林公式十分好用,而且,即使在 n 很小的时候,斯特林公式的取值已经十分准确。 斯特林矩阵(Stirling Matrix),该函数构造并返回斯特林矩阵,该矩阵是包含第二类斯特林数的下三角矩阵。 矩阵定义:

> stirling.matrix(4)
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    0    0    0    0
[2,]    0    1    0    0    0
[3,]    0    1    1    0    0
[4,]    0    1    3    1    0
[5,]    0    1    7    6    1

5.5 Pascal matrix

帕斯卡矩阵:由杨辉三角形表组成的矩阵称为帕斯卡(Pascal)矩阵。此函数返回 n 乘以 Pascal 矩阵。在数学中,尤其是矩阵理论和组合学,Pascal 矩阵是一个下三角矩阵,行中有二项式系数。通过对相同顺序的对称 Pascal 矩阵执行 LU 分解并返回下三角矩阵,可以容易地获得它。 帕斯卡的三角形是由数字行组成的三角形。第一行具有条目1.每个后续行通过添加前一行的相邻条目而形成,替换为 0,其中不存在相邻条目。pascal 函数通过选择与指定矩阵维度相对应的 Pascal 三角形部分来形成 Pascal 矩阵。 矩阵定义:

> pascal.matrix(4)
     [,1] [,2] [,3] [,4]
[1,]    1    0    0    0
[2,]    1    1    0    0
[3,]    1    2    1    0
[4,]    1    3    3    1

5.6 Fibonacci matrix

斐波纳契矩阵,该函数构造了从 Fibonacci 序列导出的 n + 1 平方 Fibonacci 矩阵。 计算公式:

> fibonacci.matrix(4)
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    1    0    0    0
[2,]    2    1    1    0    0
[3,]    3    2    1    1    0
[4,]    5    3    2    1    1
[5,]    8    5    3    2    1

5.7 Frobenius Matrix

Frobenius 矩阵也称为伴随矩阵,它出现在线性一阶微分方程组的解中。此函数返回一个在数值数学中有用的 Fronenius 矩阵。 矩阵定义:

> frobenius.matrix(4)
     [,1] [,2] [,3] [,4]
[1,]    0    0    0   -1
[2,]    1    0    0    4
[3,]    0    1    0   -6
[4,]    0    0    1    4

5.8 Duplication matrix

复制矩阵,当 A 是对称矩阵时,该函数构造将 vech(A)映射到 vec(A)的线性变换 D。 计算公式:

> D.matrix(3) 
      [,1] [,2] [,3] [,4] [,5] [,6]
 [1,]    1    0    0    0    0    0
 [2,]    0    1    0    0    0    0
 [3,]    0    0    1    0    0    0
 [4,]    0    1    0    0    0    0
 [5,]    0    0    0    1    0    0
 [6,]    0    0    0    0    1    0
 [7,]    0    0    1    0    0    0
 [8,]    0    0    0    0    1    0
 [9,]    0    0    0    0    0    1

5.9 K matrix

k 矩阵是由 H.matrices() 函数构造的,利用直积进行计算子列表的分量。K.matrix(r, c=r) ,返回阶数为 p=r*c 的方阵,对于 r 行 c 列的矩阵 A,计算 At(A) 的直积。 计算公式:

> K.matrix(2,3)
     [,1] [,2] [,3] [,4] [,5] [,6]
[1,]    1    0    0    0    0    0
[2,]    0    0    1    0    0    0
[3,]    0    0    0    0    1    0
[4,]    0    1    0    0    0    0
[5,]    0    0    0    1    0    0
[6,]    0    0    0    0    0    1

5.10 E Matrices

E 矩阵列表, E.matrices(n) 使得每个子列表的分量,是从 n 阶单位矩阵计算向量的外积导出的方阵。 计算公式:

> E.matrices(3)
[[1]]
[[1]][[1]]
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    0    0
[3,]    0    0    0

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

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


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

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

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


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

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

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

5.11 H Matrices

H 矩阵列表, H.matrices(r, c=r) 使得 r 阶 c 阶的子列表的分量,计算从 r 行和 c 列的单位矩阵的列向量的外积导出的方阵。

> H.matrices(2,3)
[[1]]
[[1]][[1]]
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    0    0

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

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


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

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

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

5.12 T Matrices

T 矩阵列表, T.matrices(n) 高级别列表中的组件数为 n。n 个组件中的每一个也是列表。每个子列表具有 n 个分量,每个分量是 n 阶矩阵。 计算公式:

> T.matrices(3)
[[1]]
[[1]][[1]]
     [,1] [,2] [,3]
[1,]    1    0    0
[2,]    0    0    0
[3,]    0    0    0

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

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


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

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

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


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

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

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

通过 R 语言,我们实现了对于矩阵的各种计算操作,非常方便!有了好的工具,不管是学习还是应用,都会事半功倍。本文只是列举了矩阵的操作方法,没有解释计算的推到过程,推到过程请参考线性代码的教科书。

—END—

本文分享自微信公众号 - 优雅R(elegant-r),作者:生信科技爱好者

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2020-05-10

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • 「Workshop」第十七期 奇异值分解

    奇异值分解(Singular Value Decomposition,以下简称SVD)是在机器学习领域广泛应用的算法,它不光可以用于降维算法中的特征分解,还可以...

    王诗翔呀
  • 「R」使用NMF包绘制热图

    学习文档: https://cran.r-project.org/web/packages/NMF/vignettes/heatmaps.pdf

    王诗翔呀
  • 绘图工具?你应该有‘掰弯’的能力

    写推文,我从来就不是纯粹地炫技,因为一旦有了TBtools,所有人都能轻易掌握并做得出来我推出来的所有图。

    王诗翔呀
  • 《Unity3D 实战核心技术详解》书中关于矩阵的错误

    最近一直在学习实时渲染,不免要接触线性代数。而渲染中,一定会用到矩阵,当我再次去复习我之前看的书时,发现《Unity3D 实战核心技术详解》关于矩阵就有几处错误...

    meteoric
  • 快速学习-机器学习(线性代数[矩阵])

    cwl_java
  • 彻底理解矩阵乘法

    今天的角度比较清奇,我们来讲讲矩阵的乘法。当然了,我告诉你的肯定不是大学教科书上那些填鸭式的云里雾里的计算规则,你可能将规则背下来了,但完全不理解为什么会这样。...

    米开朗基杨
  • 矩阵相乘法则和技巧

    只有当第一个矩阵(左矩阵)的列数等于第二矩阵(右矩阵)时,两矩阵才能相乘。因为得到的结果矩阵的i一行的第j个元素(Cij)是左矩阵第i行所有元素分别与右矩阵第...

    vv彭
  • 数据分析与数据挖掘 - 06线性代数

    导数是高等数学中非常重要的知识点,也是人工智能的算法应用中比较常用的一个知识,这一章我们的重点就是讲解一下导数和其求导法则。首先我们来看一下导数的基本概念:函数...

    马一特
  • 数学实验(预习)

    也可以用初等变换求逆矩阵,构造一个n行2n列的矩阵(A E),并进行初等变换,A编程单位矩阵的时候,E就变成了A的逆矩阵.

    云深无际
  • 吹弹牛皮之Unity 引擎基础 - 矩阵(三)

    上图中展示了p,q两个基向量(单位向量)绕原点旋转后得到的新基向量p'和q'。根据勾股定理有:

    用户7698595

扫码关注云+社区

领取腾讯云代金券