【R语言系列】R语言中矩阵运算

来源：互联网发布：手机淘宝评价管理网址编辑：程序博客网时间：2024/05/22 07:07

R语言中矩阵运算

1_矩阵的生成

2_矩阵的四则运算

3_矩阵的矩阵运算

4_矩阵的分解

1_1将向量定义成数组

向量只有定义了维数向量(dim属性)后才能被看作是数组.比如:

> z=1:12;

> dim(z)=c(3,4);

> z;

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

[1,] 1 4 7 10

[2,] 2 5 8 11

[3,] 3 6 9 12

注意：生成矩阵是按列排列的。

1_2用array ( )函数构造多维数组

用法为：array(data=NA,dim=length(data),dimnames=NULL)

参数描述：data：是一个向量数据。

dim：是数组各维的长度，缺省时为原向量的长度。

dimname：是数组维的名字，缺省时为空。

例子：

> x=array(1:20,dim=c(4,5))

> x

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

1_3用matrix()函数构造矩阵

函数matrix)是构造矩阵(二维数组)的函数，其构造形式为

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

其中data是一个向量数据，nro、是矩阵的行数，ncol是矩阵的列数.当byrow=TRUE时，生成矩阵的数据按行放置，缺省时相当于byrow=FALSE，数据按列放置.dimname。是数组维的名字，缺省时为空.

如构造一个3x5阶的矩阵

> A=matrix(1:15,nrow=3,byrow=TRUE)

> A

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

[1,] 1 2 3 4 5

[2,] 6 7 8 9 10

[3,] 11 12 13 14 15

2_矩阵的四则运算

可以对数组之间进行四则运算(+、一、*、/)，这时进行的是数组对应元素的四则运算。一般情况下参加运算的矩阵或者数组的维数是相同的，但也可以计算不同维的，这是要将对应的元素补足。

3_1 转置运算

对于矩阵A，函数t(A)表示矩阵A的转置，如：

> A=matrix(1:6,nrow=2);

> A;

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

[1,] 1 3 5

[2,] 2 4 6

> t(A);

[,1] [,2]

[1,] 1 2

[2,] 3 4

[3,] 5 6

3_2 求方阵的行列式

函数det()是求矩阵行列式的值，如

> det(matrix(1:4,ncol=2));

[1] -2

3_3 向量的内积

对于n维向量x，可以看成nxl阶矩阵或lxn阶矩阵。若x与y是相同

维数的向量，则x%*%Y表示x与y作内积.例如，

>x=1:5; Y=2*1:5

>x%*%y

[，1]

[1，]110

函数crossprod()是内积运算函数(表示交叉乘积)，crossprod(x,y)计算向量x与y的内积，即t(x) %*% y'。crossprod(x)表示x与x的内积.

类似地，tcrossprod(x,y)表示’x%*%t(Y)’，即x与y的外积，也称为叉积。tcrossprod(x)表示x与x作外积.如：

> x=1:5; y=2*1:5;

> crossprod(x);

[,1]

[1,] 55

> crossprod(x,y);

[,1]

[1,] 110

> tcrossprod(x);

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

[1,] 1 2 3 4 5

[2,] 2 4 6 8 10

[3,] 3 6 9 12 15

[4,] 4 8 12 16 20

[5,] 5 10 15 20 25

> tcrossprod(x,y);

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

[1,] 2 4 6 8 10

[2,] 4 8 12 16 20

[3,] 6 12 18 24 30

[4,] 8 16 24 32 40

[5,] 10 20 30 40 50

3_4 向量的外积(叉积)

设x和y是n维向量，则x%o%y表示x与y作外积.例如

> x%o%y;

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

[1,] 2 4 6 8 10

[2,] 4 8 12 16 20

[3,] 6 12 18 24 30

[4,] 8 16 24 32 40

[5,] 10 20 30 40 50

outer()是更为强大的外积运算函数，outer(x,y)计算向量二与y的外积，它等价于x %o%y

函数。outer()的一般调用格式为

outer(x，y，fun=”*”)

其中x, y矩阵(或向量)，fun是作外积运算函数，缺省值为乘法运算。函数outer()在绘制三维曲面时非常有用，它可生成一个x和y的网格。

3_5 矩阵的乘法

设A和B为两个矩阵，通常意义下的矩阵乘法是通过A%*%B来完成，crossprod(A,B)表示的是

t(A)%*%B，而tcrossprod(A,B)表示的是A%*%t(B)。最后我们通过运算知道x%*%A%*%x为二次型。

例子：

> A=array(1:9,dim=(c(3,3)))

> B=array(9:1,dim=(c(3,3)))

> A%*%B;

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

[1,] 90 54 18

[2,] 114 69 24

[3,] 138 84 30

> crossprod(A,B)==t(A)%*%B;

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

[1,] TRUE TRUE TRUE

[2,] TRUE TRUE TRUE

[3,] TRUE TRUE TRUE

> tcrossprod(A,B)==A%*%t(B);

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

[1,] TRUE TRUE TRUE

[2,] TRUE TRUE TRUE

[3,] TRUE TRUE TRUE

3_6 生成对角阵和矩阵取对角运算

函数diag()依赖于它的变量，当v是一个向量时，diag(v)表示以v的元素为对角线元素的对角阵.当M是一个矩阵时，则diag(M)表示的是取M对角线上的元素的向量.如

> v=c(1,4,5);

> diag(v);

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

[1,] 1 0 0

[2,] 0 4 0

[3,] 0 0 5

> M=array(1:9,dim=c(3,3));

> diag(M);

[1] 1 5 9

3_7 解线性方程组和求矩阵的逆矩阵

若求解线性方程组Ax=b，其命令形式为solve(A,b)，求矩阵A的逆，其命令形式为solve(A).设矩阵A=t(array(c(1:8,10),dim=c(3,3))),b<-c(1,1,1),则解方程组Ax=b的解x和求矩阵A的逆矩阵的命令如下:

> A=t(array(c(1:8,10),dim=c(3,3)));

> b=c(1,1,1);

> x=solve(A,b);

> x;

[1] -1.000000e+00 1.000000e+00 3.806634e-16

> solve(A);

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

[1,] -0.6666667 -1.333333 1

[2,] -0.6666667 3.666667 -2

[3,] 1.0000000 -2.000000 1

3_8 求矩阵的特征值与特征向量

函数eigen(Sm)是求对称矩阵Sm的特征值与特征向量，其命令形式为:ev=eigen(Sm),则ev存放着对称矩阵Sm特征值和特征向量，是由列表形式给出的，其中ev$values是Sm的特征值构成的向量，ev$vectors是Sm的特征向量构成的矩阵.如

> Sm=crossprod(A,A);

> ev=eigen(Sm);

> ev;

$values

[1] 303.19533618 0.76590739 0.03875643

$vectors

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

[1,] -0.4646675 0.833286355 0.2995295

[2,] -0.5537546 -0.009499485 -0.8326258

[3,] -0.6909703 -0.552759994 0.4658502

4_1 特征值分解

(1).定义：

对N阶方阵A，x为标量，v是非零的N维列向量，且满足Ax=xv ，则称x为矩阵A的特征值，v 是相对应于x 的特征向量。特征值的全体成为A的谱。

(2).在r中的实现：在r中利用函数eigen(A)来求矩阵的特征值和特征向量，具体的调用格式为：

以矩阵A为例说明此问题

> A=array(c(1,1,1,4,2,1,9,3,1),dim=c(3,3));

> D=eigen(A);

> D;

$values

[1] 5.8284271 -2.0000000 0.1715729

$vectors

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

[1,] -0.8597736 -9.486833e-01 0.5384820

[2,] -0.4346498 6.474883e-17 -0.7872938

[3,] -0.2680839 3.162278e-01 0.3003425

(3).特征值分解的性质：我们知道当所求的的特征向量构成的矩阵可逆时会满足solve(vectors)%*%A%*%vectors=diag(values),下面进行验证。

> solve(vectors)%*%A%*%vectors;

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

[1,] 5.828427e+00 8.339683e-16 -1.285213e-15

[2,] 1.211325e-15 -2.000000e+00 2.704000e-16

[3,] -3.471971e-16 -1.607126e-16 1.715729e-01

结果的精度还是比较高的。

4_2 矩阵的奇异值分解

函数svd(A)是对矩阵A作奇异值分解，即A =U%*%D%*%t(V)，其中U, V是正交阵，D为对角阵，也就是矩阵A的奇异值.svd(A)的返回值也是列表，svd(A)$d表示矩阵A的奇异值，即矩阵D的对角线上的元素.svd(A)$u对应的是正交阵U, svd(A) $v对应的是正交阵V.例如，

> A<-t(array(c(1:8,10),dim=c(3,3)))

> SVD=svd(A);

> SVD;

[1] 17.4125052 0.8751614 0.1968665

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

[1,] -0.2093373 0.96438514 0.1616762

[2,] -0.5038485 0.03532145 -0.8630696

[3,] -0.8380421 -0.26213299 0.4785099

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

[1,] -0.4646675 -0.833286355 0.2995295

[2,] -0.5537546 0.009499485 -0.8326258

[3,] -0.6909703 0.552759994 0.4658502

> attach(SVD);

The following object(s) are masked from 'SVD (position 3)':

d, u, v

> u%*%diag(d)%*%t(v);

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

[1,] 1 2 3

[2,] 4 5 6

[3,] 7 8 10

> A;

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

[1,] 1 2 3

[2,] 4 5 6

[3,] 7 8 10

4_3 qr分解

设A为m*n矩阵，如果存在m*m酉矩阵Q（即Q(H)Q=QQ(H)=I）和m*n阶梯形矩阵R，使得A=QR，那么此分解称为QR分解。QR分解在解决最小二乘问题、特征值计算等方面有着十分重要的作用。

#建立矩阵

> A=(array(c(1:12),dim=c(4,3)));

> A;

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

[1,] 1 5 9

[2,] 2 6 10

[3,] 3 7 11

[4,] 4 8 12

#进行矩阵分解

> QR=qr(A);QR

$qr

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

[1,] -5.4772256 -12.7801930 -2.008316e+01

[2,] 0.3651484 -3.2659863 -6.531973e+00

[3,] 0.5477226 -0.3781696 7.880925e-16

[4,] 0.7302967 -0.9124744 9.277920e-01

$rank

[1] 2

$qraux

[1] 1.182574 1.156135 1.373098

$pivot

[1] 1 2 3

attr(,"class")

[1] "qr"

#提取Q,R并验证分解的正确性。

> Q=qr.Q(QR);

> R=qr.R(QR);

> Q%*%R;

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

[1,] 1 5 9

[2,] 2 6 10

[3,] 3 7 11

[4,] 4 8 12

4_4 Schur分解

引言：

从特征值的分解中可以看出，特征值的分解是有条件的，如果特征向量不是线性无关的，那么对于一个矩阵来说便不能采用特征值分解的方法对矩阵进行分解。例如对于矩阵A=t(array(c(6,12,19,-9,-20,-33,4,9,15),dim=c(3,3))进行特征值分解有：

> A=t(array(c(6,12,19,-9,-20,-33,4,9,15),dim=c(3,3)));

> A;

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

[1,] 6 12 19

[2,] -9 -20 -33

[3,] 4 9 15

> det(A);

[1] -1

> W=eigen(A);

> W;

$values

[1] 1+0i 1-0i -1+0i

$vectors

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

[1,] -0.4082483-0i -0.4082483+0i -0.4740998+0i

[2,] 0.8164966+0i 0.8164966+0i 0.8127426+0i

[3,] -0.4082483+0i -0.4082483-0i -0.3386427+0i

> attach(W);

The following object(s) are masked from 'W (position 3)':

values, vectors

> det(vectors);

错误于determinant.matrix(x, logarithm = TRUE, ...) :

目前还不能算复数矩阵的行列式

> det(Re(vectors));

[1] -7.599489e-19

> solve(vectors)

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

[1,] 0.000000+78209959i 0.00000+78209959i -9.26965+78209959i

[2,] 0.000000-78209959i 0.00000-78209959i -9.10153-78209959i

[3,] 3.691206+ 0i 11.07362+ 0i 18.45603+ 0i

很明显vectors不是一个可逆矩阵此时进行特征值分辨这种方法便不可行，对于这种情况我们可以作Schur分解。

描述：

对于任意的方针A，其Schur分解的形式为：A=USU(H)，其中U是标准的正交矩阵（即满足UU(H)=I），S为上三角矩阵，并且对角线上的元素为A的特征值。由于此函数在包Matrix中，所以使用之前必须调入。并且注意matrix和Matrix的区别。

例子：

> A=Matrix(c(6,12,19,-9,-20,-33,4,9,15),ncol=3,byrow=TRUE);

> A;

3 x 3 Matrix of class "dgeMatrix"

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

[1,] 6 12 19

[2,] -9 -20 -33

[3,] 4 9 15

> library(Matrix);

> Sch=Schur(A, vectors=TRUE);

> Q=Sch@Q;

> Q=as.matrix(Q)

> attach(Sch);

错误于attach(Sch) : 'attach'只适用于串列，数据框和环境

> Q%*%T%*%t(Q)

3 x 3 Matrix of class "dgeMatrix"

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

[1,] 6 12 19

[2,] -9 -20 -33

[3,] 4 9 15

4_5 Cholesky分解（柯利分解）

描述：

正定矩阵：设A是n阶实系数矩阵，如果对任何非零向量 X=(x1,...xn) 都有

t(X)AX>0，就称A正定(Positive Definite)。正定矩阵在相合变换下可化为标准型，即单位矩阵。

Cholesky分解：

对任意的正定矩阵A，存在上三角矩阵R，使A=t(R)%*%R，则称为A的Cholesky分解（柯利分解）。

例子：

> #输入矩阵

> m=matrix(c(5,1,1,3),ncol=2 );

> m;

[,1] [,2]

[1,] 5 1

[2,] 1 3

> #矩阵分解

> CH=chol(m);

> #验证结果

> t(CH)%*%CH;

[,1] [,2]

[1,] 5 1

[2,] 1 3

0 0