提供了编程的基础技术教程

zl程序教程

您现在的位置是：首页 >  其他

当前栏目

R语言——矩阵运算

语言 矩阵 运算
2023-09-14 08:57:57 时间

R语言的矩阵运算

创建矩阵向量；矩阵加减，乘积；矩阵的逆；行列式的值；特征值与特征向量；QR分解；奇异值分解；广义逆；backsolve与fowardsolve函数；取矩阵的上下三角元素；向量化算子等。
1、创建向量
 x=c(1,2,3,4)

[1] 1 2 3 4
2、创建矩阵

在R中可以用函数matrix()来创建一个矩阵。 

 args(matrix)

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

 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

 matrix(1:12,nrow=4,ncol=3,byrow=T)

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

[1,] 1 2 3

[2,] 4 5 6

[3,] 7 8 9

[4,] 10 11 12
3、矩阵转置

A为m×n矩阵，求A’的转置矩阵在R中可用函数t()，例如：

 A=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

 t(A)

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

[1,] 1 2 3

[2,] 4 5 6

[3,] 7 8 9

[4,] 10 11 12

[1] 1 2 3 4

 x=c(1,2,3,4,5,6,7,8,9,10)

 [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] "numeric"

 class(t(x))

[1] "matrix"

 t(t(x))

 [,1]

 [1,] 1

 [2,] 2

 [3,] 3

 [4,] 4

 [5,] 5

 [6,] 6

 [7,] 7

 [8,] 8

 [9,] 9

[10,] 10

 y=t(t(x))

 t(t(y))

 [,1]

 [1,] 1

 [2,] 2

 [3,] 3

 [4,] 4

 [5,] 5

 [6,] 6

 [7,] 7

 [8,] 8

 [9,] 9

[10,] 10
4、矩阵相加减
 A=B=matrix(1:12,nrow=3,ncol=4)

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

[1,] 2 8 14 20

[2,] 4 10 16 22

[3,] 6 12 18 24

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

 [,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可用符号：“％*％”，例如：

 c=2

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

[1,] 2 8 14 20

[2,] 4 10 16 22

[3,] 6 12 18 24

 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)

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

[1,] 1 5 9 13

[2,] 2 6 10 14

[3,] 3 7 11 15

[4,] 4 8 12 16

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

[1,] 1 25 81 169

[2,] 4 36 100 196

[3,] 9 49 121 225

[4,] 16 64 144 256
7、矩阵对角元素相关运算

例如要取一个方阵的对角元素，

 A=matrix(1:16,nrow=4,ncol=4)

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

 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)

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

[1,] 0.71928674 0.4029735 0.3695724 -0.8464934

[2,] -1.06569049 0.4087710 0.8507104 0.5379580

[3,] 0.06346143 0.5549962 1.5030082 -1.2253291

[4,] 1.60231999 0.5628075 1.3339055 -1.6211637

 solve(a)

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

[1,] -0.3641840 0.2762240 -0.96264575 1.0094194

[2,] 3.4975449 1.2420380 -0.93560875 -0.7069340

[3,] -1.7608293 0.3153284 0.03335861 0.9988433

[4,] -0.5945571 0.9636571 -1.24881708 0.9572800

 solve(a)%*%a

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

[1,] 1.000000e+00 -1.110223e-16 0.000000e+00 0.000000e+00

[2,] -6.661338e-16 1.000000e+00 -6.661338e-16 8.881784e-16

[3,] 0.000000e+00 2.220446e-16 1.000000e+00 -4.440892e-16

[4,] 0.000000e+00 -1.110223e-16 2.220446e-16 1.000000e+00
9、矩阵的特征值和特征向量

矩阵A的谱分解为A=UΛU’,其中Λ是由A的特征值组成的对角矩阵，U的列为A的特征值对应的特征向量，在R中可以用函数eigen()函数得到U和Λ，

 args(eigen)

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

 A=diag(4)+1

 [,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.0000000 0.0000000

[2,] -0.5 -0.2886751 -0.5773503 -0.5773503

[3,] -0.5 -0.2886751 -0.2113249 0.7886751

[4,] -0.5 -0.2886751 0.7886751 -0.2113249

 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 -5.551115e-17 -1.110223e-16 -9.714451e-17

[2,] -5.551115e-17 1.000000e+00 -5.551115e-17 -5.551115e-17

[3,] -1.110223e-16 -5.551115e-17 1.000000e+00 0.000000e+00

[4,] -9.714451e-17 -5.551115e-17 0.000000e+00 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

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

[1] 5

 det(A)

[1] 5

 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)

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

[1] 4.589453e+01 1.640705e+00 1.366522e-15

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

[1,] -0.5290354 0.74394551 0.4082483

[2,] -0.5760715 0.03840487 -0.8164966

[3,] -0.6231077 -0.66713577 0.4082483

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

[1,] -0.07736219 -0.71960032 -0.4076688

[2,] -0.19033085 -0.50893247 0.5745647

[3,] -0.30329950 -0.29826463 -0.0280114

[4,] -0.41626816 -0.08759679 0.2226621

[5,] -0.52923682 0.12307105 -0.6212052

[6,] -0.64220548 0.33373889 0.2596585

 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 3.330669e-16 1.665335e-16

[2,] 3.330669e-16 1.000000e+00 5.551115e-17

[3,] 1.665335e-16 5.551115e-17 1.000000e+00

 t(A.svd$v)%*%A.svd$v

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

[1,] 1.000000e+00 2.775558e-17 2.775558e-17

[2,] 2.775558e-17 1.000000e+00 -2.081668e-16

[3,] 2.775558e-17 -2.081668e-16 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)

 [,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 1.601186e-15 2.217027e-15

[4,] 0.7302967 -0.9124744 -5.547002e-01 -1.478018e-15

$rank

[1] 2

$qraux

[1] 1.182574e+00 1.156135e+00 1.832050e+00 1.478018e-15

$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 1.601186e-15 2.217027e-15

[4,] 0.000000 0.000000 0.000000e+00 -1.478018e-15

 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 -1.665335e-16 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 -5.551115e-17 0.000000e+00 2.081668e-17

[2,] -5.551115e-17 1.000000e+00 -2.775558e-17 -6.938894e-17

[3,] 0.000000e+00 -2.775558e-17 1.000000e+00 2.775558e-17

[4,] 2.081668e-17 -6.938894e-17 2.775558e-17 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、矩阵的广义逆

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逆，例如：

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)

 [,1] [,2]

[1,] 1 3

[2,] 2 4

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

 [,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
16 矩阵的行和、列和、行平均与列平均

在R中很容易求得一个矩阵的各行的和、平均数与列的和、平均数，例如：

 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

上述关于矩阵行和列的操作，还可以使用apply()函数实现。
args(apply)
function (X, MARGIN, FUN, …)
其中：x为矩阵，MARGIN用来指定是对行运算还是对列运算，MARGIN＝1表示对行运算，MARGIN＝2表示对列运算，FUN用来指定运算函数, …用来给定FUN中需要的其它的参数，例如：
apply(A,1,sum)
[1] 22 26 30
apply(A,1,mean)
[1] 5.5 6.5 7.5
apply(A,2,sum)
[1] 6 15 24 33
apply(A,2,mean)
[1] 2 5 8 11
apply()函数功能强大，我们可以对矩阵的行或者列进行其它运算，例如：
计算每一列的方差
A=matrix(rnorm(100),20,5)
apply(A,2,var)
[1] 0.4641787 1.4331070 0.3186012 1.3042711 0.5238485
apply(A,2,function(x,a)x*a,a=2)
[,1] [,2] [,3] [,4]
[1,] 2 8 14 20
[2,] 4 10 16 22
[3,] 6 12 18 24
注意：apply(A,2,function(x,a)x*a,a=2)与A*2效果相同，此处旨在说明如何应用alpply函数。

17 矩阵X’X的逆

在统计计算中，我们常常需要计算这样矩阵的逆，如OLS估计中求系数矩阵。R中的包“strucchange”提供了有效的计算方法。
args(solveCrossprod)
function (X, method = c(“qr”, “chol”, “solve”))
其中：method指定求逆方法，选用“qr”效率最高，选用“chol”精度最高，选用“slove”与slove(crossprod(x,x))效果相同，例如：
A=matrix(rnorm(16),4,4)
solveCrossprod(A,method=”qr”)
[,1] [,2] [,3] [,4]
[1,] 0.6132102 -0.1543924 -0.2900796 0.2054730
[2,] -0.1543924 0.4779277 0.1859490 -0.2097302
[3,] -0.2900796 0.1859490 0.6931232 -0.3162961
[4,] 0.2054730 -0.2097302 -0.3162961 0.3447627
solveCrossprod(A,method=”chol”)
[,1] [,2] [,3] [,4]
[1,] 0.6132102 -0.1543924 -0.2900796 0.2054730
[2,] -0.1543924 0.4779277 0.1859490 -0.2097302
[3,] -0.2900796 0.1859490 0.6931232 -0.3162961
[4,] 0.2054730 -0.2097302 -0.3162961 0.3447627
solveCrossprod(A,method=”solve”)
[,1] [,2] [,3] [,4]
[1,] 0.6132102 -0.1543924 -0.2900796 0.2054730
[2,] -0.1543924 0.4779277 0.1859490 -0.2097302
[3,] -0.2900796 0.1859490 0.6931232 -0.3162961
[4,] 0.2054730 -0.2097302 -0.3162961 0.3447627
solve(crossprod(A,A))
[,1] [,2] [,3] [,4]
[1,] 0.6132102 -0.1543924 -0.2900796 0.2054730
[2,] -0.1543924 0.4779277 0.1859490 -0.2097302
[3,] -0.2900796 0.1859490 0.6931232 -0.3162961
[4,] 0.2054730 -0.2097302 -0.3162961 0.3447627

18 取矩阵的上、下三角部分

在R中，我们可以很方便的取到一个矩阵的上、下三角部分的元素，函数lower.tri()和函数upper.tri()提供了有效的方法。
args(lower.tri)
function (x, diag = FALSE)
函数将返回一个逻辑值矩阵，其中下三角部分为真，上三角部分为假，选项diag为真时包含对角元素，为假时不包含对角元素。upper.tri()的效果与之孑然相反。例如：
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
lower.tri(A)
[,1] [,2] [,3] [,4]
[1,] FALSE FALSE FALSE FALSE
[2,] TRUE FALSE FALSE FALSE
[3,] TRUE TRUE FALSE FALSE
[4,] TRUE TRUE TRUE FALSE
lower.tri(A,diag=T)
[,1] [,2] [,3] [,4]
[1,] TRUE FALSE FALSE FALSE
[2,] TRUE TRUE FALSE FALSE
[3,] TRUE TRUE TRUE FALSE
[4,] TRUE TRUE TRUE TRUE
upper.tri(A)
[,1] [,2] [,3] [,4]
[1,] FALSE TRUE TRUE TRUE
[2,] FALSE FALSE TRUE TRUE
[3,] FALSE FALSE FALSE TRUE
[4,] FALSE FALSE FALSE FALSE
upper.tri(A,diag=T)
[,1] [,2] [,3] [,4]
[1,] TRUE TRUE TRUE TRUE
[2,] FALSE TRUE TRUE TRUE
[3,] FALSE FALSE TRUE TRUE
[4,] FALSE FALSE FALSE TRUE
A[lower.tri(A)]=0
A
[,1] [,2] [,3] [,4]
[1,] 1 5 9 13
[2,] 0 6 10 14
[3,] 0 0 11 15
[4,] 0 0 0 16
A[upper.tri(A)]=0
A
[,1] [,2] [,3] [,4]
[1,] 1 0 0 0
[2,] 2 6 0 0
[3,] 3 7 11 0
[4,] 4 8 12 16

19 backsolve fowardsolve函数
这两个函数用于解特殊线性方程组，其特殊之处在于系数矩阵为上或下三角。

 args(backsolve)

function (r, x, k = ncol(r), upper.tri = TRUE, transpose = FALSE)

 args(forwardsolve)

function (l, x, k = ncol(l), upper.tri = FALSE, transpose = FALSE)

其中：r或者l为n×n维三角矩阵，x为n×1维向量，对给定不同的upper.tri和transpose的值，方程的形式不同

对于函数backsolve()而言，

 A=matrix(1:9,3,3)

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

[1,] 1 4 7

[2,] 2 5 8

[3,] 3 6 9

 x=c(1,2,3)

[1] 1 2 3

 B[upper.tri(B)]=0

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

[1,] 1 0 0

[2,] 2 5 0

[3,] 3 6 9

 C[lower.tri(C)]=0

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

[1,] 1 4 7

[2,] 0 5 8

[3,] 0 0 9

 backsolve(A,x,upper.tri=T,transpose=T)

[1] 1.00000000 -0.40000000 -0.08888889

 solve(t(C),x)

[1] 1.00000000 -0.40000000 -0.08888889

 backsolve(A,x,upper.tri=T,transpose=F)

[1] -0.8000000 -0.1333333 0.3333333

 solve(C,x)

[1] -0.8000000 -0.1333333 0.3333333

 backsolve(A,x,upper.tri=F,transpose=T)

[1] 1.111307e-17 2.220446e-17 3.333333e-01

 solve(t(B),x)

[1] 1.110223e-17 2.220446e-17 3.333333e-01

 backsolve(A,x,upper.tri=F,transpose=F)

[1] 1 0 0

 solve(B,x)

[1] 1.000000e+00 -1.540744e-33 -1.850372e-17

对于函数forwardsolve()而言，

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

[1,] 1 4 7

[2,] 2 5 8

[3,] 3 6 9

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

[1,] 1 0 0

[2,] 2 5 0

[3,] 3 6 9

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

[1,] 1 4 7

[2,] 0 5 8

[3,] 0 0 9

[1] 1 2 3

 forwardsolve(A,x,upper.tri=T,transpose=T)

[1] 1.00000000 -0.40000000 -0.08888889

 solve(t(C),x)

[1] 1.00000000 -0.40000000 -0.08888889

 forwardsolve(A,x,upper.tri=T,transpose=F)

[1] -0.8000000 -0.1333333 0.3333333

 solve(C,x)

[1] -0.8000000 -0.1333333 0.3333333

 forwardsolve(A,x,upper.tri=F,transpose=T)

[1] 1.111307e-17 2.220446e-17 3.333333e-01

 solve(t(B),x)

[1] 1.110223e-17 2.220446e-17 3.333333e-01

 forwardsolve(A,x,upper.tri=F,transpose=F)

[1] 1 0 0

 solve(B,x)

[1] 1.000000e+00 -1.540744e-33 -1.850372e-17
20 row()与col()函数
在R中定义了的这两个函数用于取矩阵元素的行或列下标矩阵，例如矩阵A={aij}m×n，

row()函数将返回一个与矩阵A有相同维数的矩阵，该矩阵的第i行第j列元素为i，函数col()类似。例如：

 x=matrix(1:12,3,4)

 row(x)

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

[1,] 1 1 1 1

[2,] 2 2 2 2

[3,] 3 3 3 3

 col(x)

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

[1,] 1 2 3 4

[2,] 1 2 3 4

[3,] 1 2 3 4

这两个函数同样可以用于取一个矩阵的上下三角矩阵，例如：

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

[1,] 1 4 7 10

[2,] 2 5 8 11

[3,] 3 6 9 12

 x[row(x) col(x)]=0

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

[1,] 1 0 0 0

[2,] 2 5 0 0

[3,] 3 6 9 0

 x=matrix(1:12,3,4)

 x[row(x) col(x)]=0

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

[1,] 1 4 7 10

[2,] 0 5 8 11

[3,] 0 0 9 12
21 行列式的值
在R中，函数det(x)将计算方阵x的行列式的值，例如：

 x=matrix(rnorm(16),4,4)

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

[1,] -1.0736375 0.2809563 -1.5796854 0.51810378

[2,] -1.6229898 -0.4175977 1.2038194 -0.06394986

[3,] -0.3989073 -0.8368334 -0.6374909 -0.23657088

[4,] 1.9413061 0.8338065 -1.5877162 -1.30568465

 det(x)

[1] 5.717667
22 向量化算子
在R中可以很容易的实现向量化算子，例如：

vec -function (x){

 t(t(as.vector(x)))

vech -function (x){

 t(x[lower.tri(x,diag=T)])

 x=matrix(1:12,3,4)

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

[1,] 1 4 7 10

[2,] 2 5 8 11

[3,] 3 6 9 12

 vec(x)

 [,1]

[1,] 1

[2,] 2

[3,] 3

[4,] 4

[5,] 5

[6,] 6

[7,] 7

[8,] 8

[9,] 9

[10,] 10

[11,] 11

[12,] 12

 vech(x)

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

[1,] 1 2 3 5 6 9
23 时间序列的滞后值
在时间序列分析中，我们常常要用到一个序列的滞后序列，R中的包“fMultivar”中的函数tslag()提供了这个功能。

 args(tslag)

function (x, k = 1, trim = FALSE)

其中：x为一个向量，k指定滞后阶数，可以是一个自然数列，若trim为假，则返回序列与原序列长度相同，但含有NA值；若trim项为真，则返回序列中不含有NA值，例如：

 x=1:20

 tslag(x,1:4,trim=F)

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

[1,] NA NA NA NA

[2,] 1 NA NA NA

[3,] 2 1 NA NA

[4,] 3 2 1 NA

[5,] 4 3 2 1

[6,] 5 4 3 2

[7,] 6 5 4 3

[8,] 7 6 5 4

[9,] 8 7 6 5

[10,] 9 8 7 6

[11,] 10 9 8 7

[12,] 11 10 9 8

[13,] 12 11 10 9

[14,] 13 12 11 10

[15,] 14 13 12 11

[16,] 15 14 13 12

[17,] 16 15 14 13

[18,] 17 16 15 14

[19,] 18 17 16 15

[20,] 19 18 17 16

 tslag(x,1:4,trim=T)

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

[1,] 4 3 2 1

[2,] 5 4 3 2

[3,] 6 5 4 3

[4,] 7 6 5 4

[5,] 8 7 6 5

[6,] 9 8 7 6

[7,] 10 9 8 7

[8,] 11 10 9 8

[9,] 12 11 10 9

[10,] 13 12 11 10

[11,] 14 13 12 11

[12,] 15 14 13 12

[13,] 16 15 14 13

[14,] 17 16 15 14

[15,] 18 17 16 15

[16,] 19 18 17 16

《R语言编程艺术》——3.2 一般矩阵运算 本节书摘来自华章计算机《R语言编程艺术》一书中的第3章，第3.2节，作者：（美）麦特洛夫（Matloff，N．）著， 更多章节内容可以访问云栖社区“华章计算机”公众号查看。
http://bbs.sciencenet.cn/home.php?mod=space&uid=443073&do=blog&id=321347 主要包括以下内容： 创建矩阵向量；矩阵加减，乘积；矩阵的逆；行列式的值；特征值与特征向量；QR分解；奇异值分解；广义逆；backsolve与fowardsolve函数；取矩阵的上下三角元素；向量化算子等.

相关文章