R语言中与矩阵相关的所有操作(上)

主要包括以下内容:
创建矩阵向量;矩阵加减,乘积;矩阵的逆;行列式的值;特征值与特征向量;QR分解;奇异值分解;广义逆;backsolve与fowardsolve函数;取矩阵的上下三角元素;向量化算子等.

1   创建一个向量
在R中可以用函数c()来创建一个向量,例如:
> x=c(1,2,3,4)
> x
[1] 1 2 3 4
 

2  
 创建一个矩阵
在R中可以用函数matrix()来创建一个矩阵,应用该函数时需要输入必要的参数值。
> args(matrix)
function (data = NA, nrow = 1, ncol = 1, byrow = FALSE, dimnames = NULL)
 

data项为必要的矩阵元素,nrow为行数,ncol为列数,注意nrowncol的乘积应为矩阵元素个数,byrow项控制排列元素时是否按行进行,dimnames给定行和列的名称。例如:
> matrix(1:12,nrow=3,ncol=4)
 
   [,1] [,2] [,3] [,4]
[1,]   1   4   7   10
[2,]   2   5   8   11
[3,]   3   6   9   12

> rowname
[1] "r1" "r2" "r3"
> colname=c("c1","c2","c3","c4")
> colname
[1] "c1" "c2" "c3" "c4"
> matrix(1:12,nrow=3,ncol=4,dimnames=list(rowname,colname))
  c1 c2 c3 c4
r1 1 4 7 10
r2 2 5 8 11

3   矩阵转置
A为m×n矩阵,求A'在R中可用函数t(),例如:
> 
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
> t(A)
   [,1] [,2] [,3]
[1,]   1   2   3
[2,]   4   5   6
[3,]   7   8   9
[4,]   10   11   12
若将函数t()作用于一个向量x,则R默认x为列向量,返回结果为一个行向量,例如:
> 
x
[1] 1 2 3 4 5 6 7 8 9 10
> t(x)
  [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,]   1   2   3   4   5   6   7   8   9   10
> class(x)
[1] "integer"
> class(t(x))
[1] "matrix"

若想得到一个列向量,可用t(t(x)),例如:
> 
x
[1] 1 2 3 4 5 6 7 8 9 10
> t(t(x))
    [,1]
[1,]   1
[2,]   2
[3,]   3
[4,]   4
[5,]   5
[6,]   6
[7,]   7
[8,]   8
[9,]   9
[10,]  10

4   矩阵相加减
在R中对同行同列矩阵相加减,可用符号:“+”、“-”,例如:
> 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
5   数与矩阵相乘
A为m×n矩阵,c>0,在R中求cA可用符号:“*”,例如:
> c=2
> c*A
 
   [,1] [,2] [,3] [,4]
[1,]   2   8   14   20
[2,]   4   10  16   22
[3,]   6   12  18   24

6  
 矩阵相乘
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

若A为n×m矩阵,要得到A'B,可用函数crossprod(),该函数计算结果与t(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)
    [,1] [,2] [,3]
[1,]  30  70 110
[2,]  70 174 278
[3,] 110 278 446
矩阵Hadamard积:若A={aij}m×n, B={bij}m×n, 则矩阵的Hadamard积定义为:
A⊙B={aij bij }m×n,R中Hadamard积可以直接运用运算符“*”例如:
> A=matrix(1:16,4,4)
> 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
> B=A
> A*B
    [,1] [,2] [,3] [,4]
[1,]   1   25   81 169
[2,]   4   36 100 196
[3,]   9   49 121 225
[4,]   16   64 144 256
R中这两个运算符的区别区加以注意。
7   矩阵对角元素相关运算
例如要取一个方阵的对角元素,
> A=matrix(1:16,nrow=4,ncol=4)
> 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
> diag(A)
[1] 1 6 11 16
对一个向量应用diag()函数将产生以这个向量为对角元素的对角矩阵,例如:
> diag(diag(A))
 
   [,1] [,2] [,3] [,4]
[1,]   1   0   0   0
[2,]   0   6   0   0
[3,]   0   0   11   0
[4,]   0   0   0   16

对一个正整数z应用diag()函数将产生以z维单位矩阵,例如:
> diag(3)
 
   [,1] [,2] [,3]
[1,]   1   0   0
[2,]   0   1   0
[3,]   0   0   1

8  
 矩阵求逆
矩阵求逆可用函数solve(),应用solve(a, b)运算结果是解线性方程组ax = b,若b缺省,则系统默认为单位矩阵,因此可用其进行矩阵求逆,例如:
> a=matrix(rnorm(16),4,4)
> a
 
           [,1]     [,2]     [,3]     [,4]
[1,] 1.6986019   0.5239738 0.2332094 0.3174184
[2,] -0.2010667 1.0913013 -1.2093734   0.8096514
[3,] -0.1797628 -0.7573283 0.2864535 1.3679963
[4,] -0.2217916 -0.3754700 0.1696771 -1.2424030
> solve(a)
              [,1]     [,2]     [,3]     [,4]
[1,] 0.9096360 0.54057479 0.7234861 1.3813059
[2,] -0.6464172 -0.91849017 -1.7546836 -2.6957775
[3,] -0.7841661 -1.78780083 -1.5795262 -3.1046207
[4,] -0.0741260 -0.06308603 0.1854137 -0.6607851
9   矩阵的特征值与特征向量
矩阵A的谱分解为A=UΛU',其中Λ是由A的特征值组成的对角矩阵,U的列为A的特征值对应的特征向量,在R中可以用函数eigen()函数得到U和Λ,
> args(eigen)
function (x, symmetric, only.values = FALSE, EISPACK = FALSE)
其中:x为矩阵,symmetric项指定矩阵x是否为对称矩阵,若不指定,系统将自动检测x是否为对称矩阵。例如:
> A=diag(4)+1
> 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
> A.eigen=eigen(A,symmetric=T)
> A.eigen
$values
[1] 5 1 1 1

$vectors
        [,1]     [,2]       [,3]     [,4]
[1,] 0.5 0.8660254 0.000000e+00 0.0000000
[2,] 0.5 -0.2886751 -6.408849e-17 0.8164966
[3,] 0.5 -0.2886751 -7.071068e-01 -0.4082483
[4,] 0.5 -0.2886751 7.071068e-01 -0.4082483

> A.eigen$vectors%*%diag(A.eigen$values)%*%t(A.eigen$vectors)
  [,1] [,2] [,3] [,4]
[1,]   2   1   1   1
[2,]   1   2   1   1
[3,]   1   1   2   1
[4,]   1   1   1   2
> t(A.eigen$vectors)%*%A.eigen$vectors
            [,1]       [,2]         [,3]         [,4]
[1,] 1.000000e+00 4.377466e-17 1.626303e-17 -5.095750e-18
[2,] 4.377466e-17 1.000000e+00 -1.694066e-18 6.349359e-18
[3,] 1.626303e-17 -1.694066e-18 1.000000e+00 -1.088268e-16
[4,] -5.095750e-18 6.349359e-18 -1.088268e-16 1.000000e+00
10   矩阵的Choleskey分解
  对于正定矩阵A,可对其进行Choleskey分解,即:A=P'P,其中P为上三角矩阵,在R中可以用函数chol()进行Choleskey分解,例如:
> 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
若矩阵为对称正定矩阵,可以利用Choleskey分解求行列式的值,如:
> prod(diag(chol(A))^2)
[1] 5
> det(A)
[1] 5
若矩阵为对称正定矩阵,可以利用Choleskey分解求矩阵的逆,这时用函数chol2inv(),这种用法更有效。如:
> chol2inv(chol(A))
      [,1] [,2] [,3] [,4]
[1,] 0.8 -0.2 -0.2 -0.2
[2,] -0.2 0.8 -0.2 -0.2
[3,] -0.2 -0.2 0.8 -0.2
[4,] -0.2 -0.2 -0.2 0.8
> solve(A)
  [,1] [,2] [,3] [,4]
[1,] 0.8 -0.2 -0.2 -0.2
[2,] -0.2 0.8 -0.2 -0.2
[3,] -0.2 -0.2 0.8 -0.2
[4,] -0.2 -0.2 -0.2 0.8

11  
 矩阵奇异值分解
  A为m×n矩阵,rank(A)= r, 可以分解为:A=UDV',其中U'U=V'V=I。在R中可以用函数scd()进行奇异值分解,例如:
> 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 3.627301e-16
  $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.18918124
[2,] -0.19033085 -0.5089325 0.42405898
[3,] -0.30329950 -0.2982646 -0.45330031
[4,] -0.41626816 -0.0875968 -0.01637004
[5,] -0.52923682 0.1230711 0.64231130
[6,] -0.64220548 0.3337389 -0.40751869
> A.svd=svd(A)
> A.svd$u%*%diag(A.svd$d)%*%t(A.svd$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
> t(A.svd$u)%*%A.svd$u
            [,1]       [,2]       [,3]
[1,] 1.000000e+00 -1.169312e-16 -3.016793e-17
[2,] -1.169312e-16 1.000000e+00 -3.678156e-17
[3,] -3.016793e-17 -3.678156e-17 1.000000e+00
> t(A.svd$v)%*%A.svd$v
        [,1]       [,2]       [,3]
[1,] 1.000000e+00 8.248068e-17 -3.903128e-18
[2,] 8.248068e-17 1.000000e+00 -2.103352e-17
[3,] -3.903128e-18 -2.103352e-17 1.000000e+00
12   矩阵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
13   矩阵广义逆(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

14  
 矩阵Kronecker积
  n×m矩阵A与h×k矩阵B的kronecker积为一个nh×mk维矩阵,
在R中kronecker积可以用函数kronecker()来计算,例如:
> A=matrix(1:4,2,2)
> B=matrix(rep(1,4),2,2)
> A
  [,1] [,2]
[1,]   1   3
[2,]   2   4
> 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
15   矩阵的维数
  在R中很容易得到一个矩阵的维数,函数dim()将返回一个矩阵的维数,nrow()返回行数,ncol()返回列数,例如:
  > 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
> nrow(A)
[1] 3
> ncol(A)
[1] 4

——————————————

往期精彩:

640?wx_fmt=png

  • 7
    点赞
  • 99
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值