总结R语言中矩阵运算的函数.docx

总结R语言中矩阵运算的函数.docx

《总结R语言中矩阵运算的函数.docx》由会员分享，可在线阅读，更多相关《总结R语言中矩阵运算的函数.docx（22页珍藏版）》请在冰豆网上搜索。

总结R语言中矩阵运算的函数

总结R语言中矩阵运算的函数

1创建一个向量

在R中可以用函数c（）来创建一个向量，例如：

>x=c（1,2,3,4）

>x

[1]1234

2创建一个矩阵

在R中可以用函数matrix（）来创建一个矩阵，应用该函数时需要输入必要的参数值。

>args（matrix）

function（data=NA,nrow=1,ncol=1,byrow=FALSE,dimnames=NULL）

data项为必要的矩阵元素，nrow为行数，ncol为列数，注意nrow与ncol的乘积应为矩阵元素个数，byrow项控制排列元素时是否按行进行，dimnames给定行和列的名称。

例如：

>matrix（1:

12,nrow=3,ncol=4）

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

[1,]14710

[2,]25811

[3,]36912

>matrix（1:

12,nrow=4,ncol=3）

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

[1,]159

[2,]2610

[3,]3711

[4,]4812

>matrix（1:

12,nrow=4,ncol=3,byrow=T）

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

[1,]123

[2,]456

[3,]789

[4,]101112

>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））

c1c2c3c4

r114710

r225811

3矩阵转置

A为m×n矩阵，求A'在R中可用函数t（），例如：

>A=matrix（1:

12,nrow=3,ncol=4）

>A

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

[1,]14710

[2,]25811

[3,]36912

>t（A）

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

[1,]123

[2,]456

[3,]789

[4,]101112

若将函数t（）作用于一个向量x，则R默认x为列向量，返回结果为一个行向量，例如：

>x

[1]12345678910

>t（x）

[,1][,2][,3][,4][,5][,6][,7][,8][,9][,10]

[1,]12345678910

>class（x）

[1]"integer"

>class（t（x））

[1]"matrix"

若想得到一个列向量，可用t（t（x）），例如：

>x

[1]12345678910

>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矩阵相加减

在R中对同行同列矩阵相加减，可用符号：

“＋”、“－”，例如：

>A=B=matrix（1:

12,nrow=3,ncol=4）

>A+B

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

[1,]281420

[2,]4101622

[3,]6121824

>A-B

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

[1,]0000

[2,]0000

[3,]0000

5数与矩阵相乘

A为m×n矩阵，c>0，在R中求cA可用符号：

“*”，例如：

>c=2

>c*A

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

[1,]281420

[2,]4101622

[3,]6121824

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,]70158246

[2,]80184288

[3,]90210330

若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,]3070110

[2,]70174278

[3,]110278446

>crossprod（A,B）

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

[1,]3070110

[2,]70174278

[3,]110278446

矩阵Hadamard积：

若A={aij}m×n,B={bij}m×n,则矩阵的Hadamard积定义为：

A⊙B={aijbij}m×n,R中Hadamard积可以直接运用运算符“*”例如：

>A=matrix（1:

16,4,4）

>A

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

[1,]15913

[2,]261014

[3,]371115

[4,]481216

>B=A

>A*B

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

[1,]12581169

[2,]436100196

[3,]949121225

[4,]1664144256

R中这两个运算符的区别区加以注意。

7矩阵对角元素相关运算

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

>A=matrix（1:

16,nrow=4,ncol=4）

>A

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

[1,]15913

[2,]261014

[3,]371115

[4,]481216

>diag（A）

[1]161116

对一个向量应用diag（）函数将产生以这个向量为对角元素的对角矩阵，例如：

>diag（diag（A））

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

[1,]1000

[2,]0600

[3,]00110

[4,]00016

对一个正整数z应用diag（）函数将产生以z维单位矩阵，例如：

>diag（3）

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

[1,]100

[2,]010

[3,]001

8矩阵求逆

矩阵求逆可用函数solve（），应用solve（a,b）运算结果是解线性方程组ax=b，若b缺省，则系统默认为单位矩阵，因此可用其进行矩阵求逆，例如：

>a=matrix（rnorm（16）,4,4）

>a

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

[1,]1.69860190.52397380.23320940.3174184

[2,]-0.20106671.0913013-1.20937340.8096514

[3,]-0.1797628-0.75732830.28645351.3679963

[4,]-0.2217916-0.37547000.1696771-1.2424030

>solve（a）

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

[1,]0.90963600.540574790.72348611.3813059

[2,]-0.6464172-0.91849017-1.7546836-2.6957775

[3,]-0.7841661-1.78780083-1.5795262-3.1046207

[4,]-0.0741260-0.063086030.1854137-0.6607851

>solve（a）%*%a

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

[1,]1.000000e+002.748453e-17-2.787755e-17-8.023096e-17

[2,]1.626303e-191.000000e+00-4.960225e-186.977925e-16

[3,]2.135878e-17-4.629543e-171.000000e+006.201636e-17

[4,]1.866183e-171.563962e-171.183813e-171.000000e+00

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,]2111

[2,]1211

[3,]1121

[4,]1112

>A.eigen=eigen（A,symmetric=T）

>A.eigen

$values

[1]5111

$vectors

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

[1,]0.50.86602540.000000e+000.0000000

[2,]0.5-0.2886751-6.408849e-170.8164966

[3,]0.5-0.2886751-7.071068e-01-0.4082483

[4,]0.5-0.28867517.071068e-01-0.4082483

>A.eigen$vectors%*%diag（A.eigen$values）%*%t（A.eigen$vectors）

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

[1,]2111

[2,]1211

[3,]1121

[4,]1112

>t（A.eigen$vectors）%*%A.eigen$vectors

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

[1,]1.000000e+004.377466e-171.626303e-17-5.095750e-18

[2,]4.377466e-171.000000e+00-1.694066e-186.349359e-18

[3,]1.626303e-17-1.694066e-181.000000e+00-1.088268e-16

[4,]-5.095750e-186.349359e-18-1.088268e-161.000000e+00

10矩阵的Choleskey分解

对于正定矩阵A，可对其进行Choleskey分解，即：

A=P'P，其中P为上三角矩阵，在R中可以用函数chol（）进行Choleskey分解，例如：

>A

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

[1,]2111

[2,]1211

[3,]1121

[4,]1112

>chol（A）

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

[1,]1.4142140.70710680.70710680.7071068

[2,]0.0000001.22474490.40824830.4082483

[3,]0.0000000.00000001.15470050.2886751

[4,]0.0000000.00000000.00000001.1180340

>t（chol（A））%*%chol（A）

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

[1,]2111

[2,]1211

[3,]1121

[4,]1112

>crossprod（chol（A）,chol（A））

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

[1,]2111

[2,]1211

[3,]1121

[4,]1112

若矩阵为对称正定矩阵，可以利用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.20.8-0.2-0.2

[3,]-0.2-0.20.8-0.2

[4,]-0.2-0.2-0.20.8

>solve（A）

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

[1,]0.8-0.2-0.2-0.2

[2,]-0.20.8-0.2-0.2

[3,]-0.2-0.20.8-0.2

[4,]-0.2-0.2-0.20.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,]147101316

[2,]258111417

[3,]369121518

>svd（A）

$d

[1]4.589453e+011.640705e+003.627301e-16

$u

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

[1,]-0.52903540.743945510.4082483

[2,]-0.57607150.03840487-0.8164966

[3,]-0.6231077-0.667135770.4082483

$v

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

[1,]-0.07736219-0.7196003-0.18918124

[2,]-0.19033085-0.50893250.42405898

[3,]-0.30329950-0.2982646-0.45330031

[4,]-0.41626816-0.0875968-0.01637004

[5,]-0.529236820.12307110.64231130

[6,]-0.642205480.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,]147101316

[2,]258111417

[3,]369121518

>t（A.svd$u）%*%A.svd$u

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

[1,]1.000000e+00-1.169312e-16-3.016793e-17

[2,]-1.169312e-161.000000e+00-3.678156e-17

[3,]-3.016793e-17-3.678156e-171.000000e+00

>t（A.svd$v）%*%A.svd$v

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

[1,]1.000000e+008.248068e-17-3.903128e-18

[2,]8.248068e-171.000000e+00-2.103352e-17

[3,]-3.903128e-18-2.103352e-171.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.37816962.641083e-152.056562e-15

[4,]0.7302967-0.91247448.583032e-01-2.111449e-16

$rank

[1]2

$qraux

[1]1.182574e+001.156135e+001.513143e+002.111449e-16

$pivot

[1]1234

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.0000000.0000002.641083e-152.056562e-15

[4,]0.0000000.0000000.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-010.25463290.79697056

[3,]-0.5477226-8.131516e-190.6909965-0.47172438

[4,]-0.73029674.082483e-01-0.54554190.04882607

>qr.Q（qr（A））%*%qr.R（qr（A））

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

[1,]15913

[2,]261014

[3,]371115

[4,]481216

>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-161.000000e+00-4.269046e-177.011739e-17

[3,]-6.760001e-17-4.269046e-171.000000e+00-1.596437e-16

[4,]-7.659550e-177.011739e-17-1.596437e-161.000000e+00

>qr.X（qr（A））

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

[1,]15913

[2,]261014

[3,]371115

[4,]481216

13矩阵广义逆（Moore-Penrose）

n×m矩阵A+称为m×n矩阵A的Moore-Penrose逆，如果它满足下列条件：

①AA+A=A；②A+AA+=A+；③（AA+）H=AA+；④（A+A）H=A+A

在R的MASS包中的函数ginv（）可计算矩阵A的Moore-Penrose逆，例如：

library（“MASS”）

>A

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

[1,]15913

[2,]261014

[3,]371115

[4,]481216

>ginv（A）

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

[1,]-0.285-0.10750.070.2475

[2,]-0.145-0.05250.040.1325

[3,]-0.0050.00250.010.0175

[4,]0.1350.0575-0.02-0.0975

验证性质1：

>A%*%ginv（A）%*%A

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

[1,]15913

[2,]261014

[3,]371115

[4,]481216

验证性质2：

>ginv（A）%*%A%*%ginv（A）

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

[1,]-0.285-0.10750.070.2475

[2,]-0.145-0.05250.040.1325

[3,]-0.0050.00250.010.0175

[4,]0.1350.0575-0.02-0.0975

验证性质3:

>t（A%*%ginv（A））

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

[1,]0.70.40.1-0.2

[2,]0.40.30.20.1

[3,]0.10.20.30.4

[4,]-0.20.10.40.7

>A%*%ginv（A）

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

[1,]0.70.40.1-0.2

[2,]0.40.30.20.1

[3,]0.10.20.30.4

[4,]-0.20.10.40.7

验证性质4:

>t（ginv（A）%*%A）

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

[1,]0.70.40.1-0.2

[2,]0.40.30.20.1

[3,]0.10.20.30.4

[4,]-0.20.10.40.7

>ginv（A）%*%A

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

[1,]0.70.40.1-0.2

[2,]0.40.30.20.1

[3,]0.10.20.30.4

[4,]-0.20.10.40.7

14矩阵Kronecker积

n×m矩阵