R中的矩阵运算-基本运算

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 为列数，注意nrow 与ncol 的乘积应为
矩阵元素个数，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
> matrix(1:12,nrow=4,ncol=3)
[,1] [,2] [,3]
[1,] 1 5 9
[2,] 2 6 10
[3,] 3 7 11
[4,] 4 8 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
> 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

也可以用array这个函数

> args(array)
function (data = NA, dim = length(data), dimnames = NULL) 

> array(c(1:6), c(2,3))
     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6

3.矩阵的维数

在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
4. 矩阵的行和、列和、行平均与列平均

在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

5. 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
这两个函数同样可以用于取一个矩阵的上下三角矩阵，例如：
> x
[,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
> x
[,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
> x
[,1] [,2] [,3] [,4]
[1,] 1 4 7 10
[2,] 0 5 8 11
[3,] 0 0 9 12

6. 取矩阵的上、下三角部分
在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

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 矩阵转置
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
> 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

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

10. 数与矩阵相乘
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

11.向量、矩阵的内积
    对于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

12.向量、矩阵的外积(叉积)
设x和y是n维向量，则x%o%y表示x与y作外积.

An important operation on arrays is the outer product. If a and b are two numeric arrays,their outer product is an array whose 

dimension vector is obtained by concatenating their two dimension vectors (order is important), and whose data vector is got 

by forming all possible products of elements of the data vector of a with those of b.

> x
[1] 1 2 3
> y
[1] 2 3

向量的维数默认是其长度；

> x %o% y
     [,1] [,2]
[1,]    2    3
[2,]    4    6
[3,]    6    9

> a
     [,1] [,2] [,3] [,4]
[1,]    1    4    7   10
[2,]    2    5    8   11
[3,]    3    6    9   12

> b
     [,1] [,2]
[1,]    1    3
[2,]    2    4
> a %o% b
, , 1, 1
     [,1] [,2] [,3] [,4]
[1,]    1    4    7   10
[2,]    2    5    8   11
[3,]    3    6    9   12

, , 2, 1
     [,1] [,2] [,3] [,4]
[1,]    2    8   14   20
[2,]    4   10   16   22
[3,]    6   12   18   24
, , 1, 2
     [,1] [,2] [,3] [,4]
[1,]    3   12   21   30
[2,]    6   15   24   33
[3,]    9   18   27   36
, , 2, 2
     [,1] [,2] [,3] [,4]
[1,]    4   16   28   40
[2,]    8   20   32   44
[3,]   12   24   36   48

      outer()是更为强大的外积运算函数，outer(x,y)计算向量x与y的外积，它等价于x %o%y
函数。outer()的一般调用格式为
      outer(x，y，fun=”*”)
     其中x, y矩阵(或向量)，fun是作外积运算函数，缺省值为乘法运算。函数outer()在绘制三维曲面时非常有用，它可生成一个x和y的网格。

13 行列式的值
在R中，函数det(x)将计算方阵x的行列式的值，例如：
> x=matrix(rnorm(16),4,4)
> x
[,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

14.矩阵的秩

> a
     [,1] [,2] [,3] [,4]
[1,]    1    4    7   10
[2,]    0    5    8   11
[3,]    0    0    9   12

> qr(a)$rank
[1] 3

15.解线性方程组和求矩阵的逆矩阵

    若求解线性方程组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

16.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
an upper (or lower) triangular matrix giving the coefficients for the system to be solved.  Values below (above) the diagonal are ignored. 
上（或下）三角矩阵给系统要解决的系数。值低于（高于）对角线被忽略。
参数：x
a matrix whose columns give the right-hand sides for the equations. 
一个矩阵的列给方程的右端。
参数：k
The number of columns of r and rows of x to use. 
r和x使用的行的列数。
参数：upper.tri
logical; if TRUE (default), the upper triangular part of r is used.  Otherwise, the lower one. 
逻辑;如果TRUE（默认），r上三角部分使用。否则，较低的一个。
参数：transpose
logical; if TRUE, solve r' * y = x for y, i.e., t(r) %*% y == x. 
如果TRUE，解决逻辑;r' * y = xy，即t(r) %*% y == x。

> A=matrix(1:9,3,3)
> A
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9
> x=c(1,2,3)
> x
[1] 1 2 3
> B=A
> B[upper.tri(B)]=0
> B
[,1] [,2] [,3]
[1,] 1 0 0
[2,] 2 5 0
[3,] 3 6 9
> C=A
> C[lower.tri(C)]=0
> C
[,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

17.特征值和特征向量

矩阵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

> eigen(a)
$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

> c <- eigen(a)

> c[[2]] %*% diag(c[[1]]) %*% solve(c[[2]])
     [,1] [,2] [,3] [,4]
[1,]    2    1    1    1
[2,]    1    2    1    1
[3,]    1    1    2    1
[4,]    1    1    1    2