1 矩阵基本操作
1.1创建向量
R里面有多种方法来创建向量（Vector），最简单的是用函数c()。例如：
>X=c(1,2,3,4)
>X
[1] 1 2 3 4
当然，还有别的方法。例如：
>X=1:4
>X
[1] 1 2 3 4
还有seq()函数。例如：
> X=seq(1,4,length=4)
> X
[1] 1 2 3 4
注意一点，R中的向量默认为列向量，如果要得到行向量需要对其进行转置。
1.2创建矩阵
R中创建矩阵的方法也有很多。大致分为直接创建和由其它格式转换两种方法。
1.2.1直接创建矩阵
最简单的直接创建矩阵的方法是用matrix()函数，matrix()函数的使用方法如下：
> args(matrix)
function (data = NA, nrow = 1, ncol = 1, byrow = FALSE, dimnames = NULL)
NULL
其中，data参数输入的为矩阵的元素，不能为空；nrow参数输入的是矩阵的行数，默认为1；ncol参数输入的是矩阵的列数，默认为1；byrow参数控制矩阵元素的排列方式，TRUE表示按行排列，FALSE表示按列排列，默认为FALSE；dimnames参数输入矩阵的行名和列名，可以不输入，系统默认为NULL。例如：
> matrix(1:6,nrow=2,ncol=3,byrow=FALSE)
[,1] [,2] [,3]
[1,] 1 3 5
[2,] 2 4 6
改变矩阵的行数和列数：
> matrix(1:6,nrow=3,ncol=2,byrow=FALSE)
[,1] [,2]
[1,] 1 4
[2,] 2 5
[3,] 3 6
改变byrow参数：
> matrix(1:6,nrow=3,ncol=2,byrow=T)
[,1] [,2]
[1,] 1 2
[2,] 3 4
[3,] 5 6
设定矩阵的行名和列名：
> matrix(1:6,nrow=3,ncol=2,byrow=T,dimnames=list(c(“A”,”B”,”C”),c(“boy”,”girl”)))
boy girl
A 1 2
B 3 4
C 5 6
1.2.2 由其它格式转换
也可以由其它格式的数据转换为矩阵，此时需要用到函数as.matrix()。
1.3 查看和改变矩阵的维数
矩阵有两个维数，即行维数和列维数。在R中查看矩阵的行维数和列维数可以用函数dim()。例如：
> X=matrix(1:12,ncol=3,nrow=4)
> X
[,1] [,2] [,3]
[1,] 1 5 9
[2,] 2 6 10
[3,] 3 7 11
[4,] 4 8 12
> dim(X)
[1] 4 3
只返回行维数：
> dim(X)[1]
[1] 4
也可以用函数nrow()
> nrow(X)
[1] 4
只返回列维数：
> dim(X)[2]
[1] 3
也可以用函数ncol():
> ncol(X)
[1] 3
同时，函数dim()也可以改变矩阵的维数。例如：
> dim(X)=c(2,6)
> X
[,1] [,2] [,3] [,4] [,5] [,6]
[1,] 1 3 5 7 9 11
[2,] 2 4 6 8 10 12
1.4矩阵行列的名称
查看矩阵的行名和列名分别用函数rownames()和函数colnames()。例如：
> X=matrix(1:6,nrow=3,ncol=2,byrow=T,dimnames=list(c(“A”,”B”,”C”),c(“boy”,”girl”)))
> X
boy girl
A 1 2
B 3 4
C 5 6
查看矩阵的行名：
> rownames(X)
[1] “A” “B” “C”
查看矩阵的列名：
> colnames(X)
[1] “boy” “girl”
同时也可以改变矩阵的行名和列名，比如：
> rownames(X)=c(“E”,”F”,”G”)
> X
boy girl
E 1 2
F 3 4
G 5 6
> colnames(X)=c(“man”,”woman”)
> X
man woman
E 1 2
F 3 4
G 5 6
1.5 矩阵元素的查看及其重新赋值
查看矩阵的某个特定元素，只需要知道该元素的行坐标和列坐标即可，例如：
> X=matrix(1:12,nrow=4,ncol=3)
> X
[,1] [,2] [,3]
[1,] 1 5 9
[2,] 2 6 10
[3,] 3 7 11
[4,] 4 8 12
查看位于矩阵第二行、第三列的元素：
> X[2,3]
[1] 10
也可以重新对矩阵的元素进行赋值，将矩阵第二行、第三列的元素替换为0：
> X[2,3]=0
> X
[,1] [,2] [,3]
[1,] 1 5 9
[2,] 2 6 0
[3,] 3 7 11
[4,] 4 8 12
R中有一个diag()函数可以返回矩阵的全部对角元素：
> X=matrix(1:9,ncol=3,nrow=3)
> X
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9
> diag(X)
[1] 1 5 9
当然也可以对对角元素进行重新赋值：
> diag(X)=c(0,0,1)
> X
[,1] [,2] [,3]
[1,] 0 4 7
[2,] 2 0 8
[3,] 3 6 1
当操作对象不是对称矩阵时，diag()也可以进行操作。
> X=matrix(1:12,nrow=4,ncol=3)
> X
[,1] [,2] [,3]
[1,] 1 5 9
[2,] 2 6 10
[3,] 3 7 11
[4,] 4 8 12
> diag(X)
[1] 1 6 11
diag()函数还能用来生成对角矩阵：
> diag(c(1,2,3))
[,1] [,2] [,3]
[1,] 1 0 0
[2,] 0 2 0
[3,] 0 0 3
也可以生成单位对角矩阵：
> diag(3)
[,1] [,2] [,3]
[1,] 1 0 0
[2,] 0 1 0
[3,] 0 0 1
> diag(4)
[,1] [,2] [,3] [,4]
[1,] 1 0 0 0
[2,] 0 1 0 0
[3,] 0 0 1 0
[4,] 0 0 0 1
查看矩阵的上三角部分：在R中查看矩阵的上三角和下三角部分很简单。可以通过lower.tri()和upper.tir()来实现：
> args(lower.tri)
function (x, diag = FALSE)
NULL
> args(upper.tri)
function (x, diag = FALSE)
NULL
查看上三角：
> X=matrix(1:12,nrow=4,ncol=3)
> X
[,1] [,2] [,3]
[1,] 1 5 9
[2,] 2 6 10
[3,] 3 7 11
[4,] 4 8 12
> X[lower.tri(X)]
[1] 2 3 4 7 8 12
改变赋值：
> X[lower.tri(X)]=0
> X
[,1] [,2] [,3]
[1,] 1 5 9
[2,] 0 6 10
[3,] 0 0 11
[4,] 0 0 0
2 矩阵计算
2.1矩阵转置
R中矩阵的转置可以用t()函数完成，例如：
> X=matrix(1:12,nrow=4,ncol=3)
> X
[,1] [,2] [,3]
[1,] 1 5 9
[2,] 2 6 10
[3,] 3 7 11
[4,] 4 8 12
> t(X)
[,1] [,2] [,3] [,4]
[1,] 1 2 3 4
[2,] 5 6 7 8
[3,] 9 10 11 12
2.2矩阵的行和与列和及行平均值和列均值
在R中很容易计算一个矩阵的各行和和各列和以及各行的平均值和各列的平均值。例如：
> 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
> 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
2.3行列式的值
R中的函数det()将计算方阵A的行列式。例如：
> X=matrix(rnorm(9),nrow=3,ncol=3)
> X
[,1] [,2] [,3]
[1,] 0.05810412 -1.2992698 0.5630315
[2,] -0.28070583 0.1958623 -1.8202283
[3,] 0.83691209 0.4411497 1.0014306
> det(X)
[1] 1.510076
2.4矩阵相加减
矩阵元素的相加减是指维数相同的矩阵，处于同行和同列的位置的元素进行加减。实现这个功能用“＋”，“－”即可。例如：
> 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
2.5矩阵的数乘
矩阵的数乘是指一个常数与一个矩阵相乘。设A为m×n矩阵，c≠0，在R中求cA的值，可以用符号“*”。例如：
> c=2
> 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
> c*A
[,1] [,2] [,3] [,4]
[1,] 2 8 14 20
[2,] 4 10 16 22
[3,] 6 12 18 24
结果矩阵与原矩阵的所有相应元素都差一个常数c。
2.6矩阵相乘
2.6.1矩阵的乘法
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
注意BA无意义，因其不符合矩阵的相乘规则。
若A为n×m矩阵，B为n×k矩阵，在R中求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：
> crossprod(A,B)
[,1] [,2] [,3]
[1,] 30 70 110
[2,] 70 174 278
[3,] 110 278 446
2.6.2矩阵的Kronecker积
n×m矩阵A和h×k矩阵B的Kronecker积是一个nh×mk维矩阵，公式为：
a11B … a1nB
Am×n×Bh×k= … …
am1B … amnB mh×nk
在R中Kronecker积可以用函数kronecher()来计算。例如：
> A=matrix(1:4,2,2)
> A
[,1] [,2]
[1,] 1 3
[2,] 2 4
> B=matrix(rep(1,4),2,2)
> 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
2.7矩阵的伴随矩阵
求矩阵A的伴随矩阵可以用LoopAnalyst包中的函数make.adjoint()函数。例如：
>install.packages(“LoopAnalyst”)
> 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
> make.adjoint(A)
[,1] [,2] [,3]
[1,] -3 6 -3
[2,] 6 -12 6
[3,] -3 6 -3
2.8矩阵的逆和广义逆
2.8.1矩阵的逆
矩阵A的逆A-1可以用函数solve()，例如：
> A=matrix(rnorm(9),nrow=3,ncol=3)
> A
[,1] [,2] [,3]
[1,] -0.2915845 0.2831544 0.94493154
[2,] -1.6494678 0.6999185 -0.06292334
[3,] -0.7224015 -0.3906971 0.44799963
> solve(A)
[,1] [,2] [,3]
[1,] 0.2359821 -0.4050650 -0.5546321
[2,] 0.6405592 0.4507583 -1.2877720
[3,] 0.9391490 -0.2600663 0.2147417
验证AA-1=1：
> A%*%solve(A)
[,1] [,2] [,3]
[1,] 1.000000e+00 8.433738e-17 -1.341700e-18
[2,] 1.216339e-17 1.000000e+00 -4.667152e-17
[3,] -2.203641e-17 4.283954e-17 1.000000e+00
用round函数可以更好的得到结果：
> round(A%*%solve(A))
[,1] [,2] [,3]
[1,] 1 0 0
[2,] 0 1 0
[3,] 0 0 1
solve()函数也可以用来求解方程组ax=b。
2.8.2矩阵的广义逆（Moore-Penrose）
并非所有的矩阵都有逆，但是所有的矩阵都可有广义逆。n×m矩阵A+是矩阵A的Moore-Penrose逆，如果它满足下列条件：
AA+A=A
A+AA+=A+
(AA+)T=AA+
(A+A)T=A+A
R中MASS包中的ginv()函数可以计算矩阵的Moore-Penrose逆。例如：
> library(MASS)
> 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
> solve(A)
Error in solve.default(A) : only square matrices can be inverted
> ginv(A)
[,1] [,2] [,3]
[1,] -0.483333333 -0.03333333 0.41666667
[2,] -0.244444444 -0.01111111 0.22222222
[3,] -0.005555556 0.01111111 0.02777778
[4,] 0.233333333 0.03333333 -0.16666667
验证性质①：
> A%*%ginv(A)%*%A
[,1] [,2] [,3] [,4]
[1,] 1 4 7 10
[2,] 2 5 8 11
[3,] 3 6 9 12
> A
[,1] [,2] [,3] [,4]
[1,] 1 4 7 10
[2,] 2 5 8 11
[3,] 3 6 9 12
验证性质②：
> ginv(A)%*%A%*%ginv(A)
[,1] [,2] [,3]
[1,] -0.483333333 -0.03333333 0.41666667
[2,] -0.244444444 -0.01111111 0.22222222
[3,] -0.005555556 0.01111111 0.02777778
[4,] 0.233333333 0.03333333 -0.16666667
> ginv(A)
[,1] [,2] [,3]
[1,] -0.483333333 -0.03333333 0.41666667
[2,] -0.244444444 -0.01111111 0.22222222
[3,] -0.005555556 0.01111111 0.02777778
[4,] 0.233333333 0.03333333 -0.16666667
验证性质③：
> A%*%ginv(A)
[,1] [,2] [,3]
[1,] 0.8333333 0.3333333 -0.1666667
[2,] 0.3333333 0.3333333 0.3333333
[3,] -0.1666667 0.3333333 0.8333333
> t(A%*%ginv(A))
[,1] [,2] [,3]
[1,] 0.8333333 0.3333333 -0.1666667
[2,] 0.3333333 0.3333333 0.3333333
[3,] -0.1666667 0.3333333 0.8333333
验证性质④：
> 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
> 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
也可以不必如此麻烦来验证性质③和④，因为③和④只是表明AA+和A+A是对称矩阵。
2.8.3 X’X的逆
很多时候，我们需要计算形如X’X的逆。这很容易实现，例如：
> x=matrix(rnorm(9),ncol=3,nrow=3)
> x
[,1] [,2] [,3]
[1,] -0.1806586 -0.76340512 0.002652331
[2,] -1.8018584 0.04467943 1.416332187
[3,] 1.2785359 -1.31653513 0.180653002
> solve(crossprod(x))
[,1] [,2] [,3]
[1,] 1.2181837 0.9664576 1.470940
[2,] 0.9664576 1.2010110 1.204599
[3,] 1.4709402 1.2045986 2.269921
R中的strucchange包中的函数solveCrossprod()也可完成：
> args(solveCrossprod)
function (X, method = c(“qr”, “chol”, “solve”))
NULL
> solveCrossprod(x,method=”qr”)
[,1] [,2] [,3]
[1,] 1.2181837 0.9664576 1.470940
[2,] 0.9664576 1.2010110 1.204599
[3,] 1.4709402 1.2045986 2.269921
> solveCrossprod(x,method=”chol”)
[,1] [,2] [,3]
[1,] 1.2181837 0.9664576 1.470940
[2,] 0.9664576 1.2010110 1.204599
[3,] 1.4709402 1.2045986 2.269921
> solveCrossprod(x,method=”solve”)
[,1] [,2] [,3]
[1,] 1.2181837 0.9664576 1.470940
[2,] 0.9664576 1.2010110 1.204599
[3,] 1.4709402 1.2045986 2.269921
2.9矩阵的特征值和特征向量
可以通过对矩阵A进行谱分解来得到矩阵的特征值和特征向量。矩阵A的谱分解如下：A=UΛU’，其中U的列为A的特征值所对应的特征向量，在R中可以用eigen()函数得到U和Λ。例如：
> args(eigen)
function (x, symmetric, only.values = FALSE, EISPACK = FALSE)
NULL
其中，x参数输入矩阵；symmetric参数判断矩阵是否为对称矩阵，如果参数为空，系统将自动检测矩阵的对称性。例如：
> A=matrix(1:9,nrow=3,ncol=3)
> A
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9
> Aeigen=eigen(A)
> Aeigen
$values
[1] 1.611684e+01 -1.116844e+00 -4.054214e-16
$vectors
[,1] [,2] [,3]
[1,] -0.4645473 -0.8829060 0.4082483
[2,] -0.5707955 -0.2395204 -0.8164966
[3,] -0.6770438 0.4038651 0.4082483
得到矩阵A的特征值：
> Aeigen$values
[1] 1.611684e+01 -1.116844e+00 -4.054214e-16
得到矩阵A的特征向量：
> Aeigen$vectors
[,1] [,2] [,3]
[1,] -0.4645473 -0.8829060 0.4082483
[2,] -0.5707955 -0.2395204 -0.8164966
[3,] -0.6770438 0.4038651 0.4082483
3 矩阵高级操作
3.1 Choleskey分解
对于正定矩阵A，可以对其进行Choleskey分解，A=P’P，其中P为上三角矩阵，在R中可以用函数chol()。例如：
> A=diag(3)+1
> A
[,1] [,2] [,3]
[1,] 2 1 1
[2,] 1 2 1
[3,] 1 1 2
> chol(A)
[,1] [,2] [,3]
[1,] 1.414214 0.7071068 0.7071068
[2,] 0.000000 1.2247449 0.4082483
[3,] 0.000000 0.0000000 1.1547005
验证A=P’P：
> t(chol(A))%*%chol(A)
[,1] [,2] [,3]
[1,] 2 1 1
[2,] 1 2 1
[3,] 1 1 2
也可以用crossprod()函数：
> crossprod(chol(A),chol(A))
[,1] [,2] [,3]
[1,] 2 1 1
[2,] 1 2 1
[3,] 1 1 2
可以用Choleskey分解来计算矩阵的行列式：
> prod(diag(chol(A))^2)
[1] 4
> det(A)
[1] 4
也可以用Choleskey分解来计算矩阵的逆，这时候可以用到函数chol2inv():
> chol2inv(chol(A))
[,1] [,2] [,3]
[1,] 0.75 -0.25 -0.25
[2,] -0.25 0.75 -0.25
[3,] -0.25 -0.25 0.75
> solve(A)
[,1] [,2] [,3]
[1,] 0.75 -0.25 -0.25
[2,] -0.25 0.75 -0.25
[3,] -0.25 -0.25 0.75
3.2奇异值分解
A为m×n矩阵，矩阵的秩为r。A可以分解为A=UDV’，其中U’U=V’V=I。在R中可以用函数svd()。例如：
> 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 2.294505e-15
$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.67039144
[2,] -0.19033085 -0.5089325 0.55766549
[3,] -0.30329950 -0.2982646 0.28189237
[4,] -0.41626816 -0.0875968 0.07320847
[5,] -0.52923682 0.1230711 0.12920119
[6,] -0.64220548 0.3337389 -0.37157608
> A.u%*%diag(A.d)%*%t(A.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
3.3 QR分解
A为m×n矩阵可以进行QR分解:A=QR，其中Q’Q=I，在R中可以用函数qr()来完成这个过程，例如：
> A=matrix(1:12,4,3)
> qr(A)
$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”
Rank返回的是矩阵的秩。Qr项包含了Q矩阵和R矩阵的信息。要想得到Q矩阵和R矩阵，可以用qr.Q()函数和qr.R()函数：
> qr.Q(qr(A))
[,1] [,2] [,3]
[1,] -0.1825742 -8.164966e-01 -0.4000874
[2,] -0.3651484 -4.082483e-01 0.2546329
[3,] -0.5477226 4.938541e-17 0.6909965
[4,] -0.7302967 4.082483e-01 -0.5455419
> qr.R(qr(A))
[,1] [,2] [,3]
[1,] -5.477226 -12.780193 -2.008316e+01
[2,] 0.000000 -3.265986 -6.531973e+00
[3,] 0.000000 0.000000 7.880925e-16
4 解方程组
4.1普通方程组
解普通方程组可以用函数solve()，solve()的基本用法是solve(A,b)，其中，A为方程组的系数矩阵，b为方程组的右端。例如：
已知方程组：
2x1+2x3=1
2x1+x2+2x3=2
2x1+x2=3
解法如下：
> A
[,1] [,2] [,3]
[1,] 2 0 2
[2,] 2 1 2
[3,] 2 1 0
> b=1:3
>b
[1] 1 2 3
> solve(A,b)
[1] 1.0 1.0 -0.5
即x1=1，x2=1，x3=-0.5。
4.2 特殊方程组
对于系数矩阵是上三角矩阵和下三角矩阵的方程组。R中提供了backsolve()和fowardsolve()来解决这个问题。
backsolve(r, x, k=ncol(r), upper.tri=TRUE, transpose=FALSE)
forwardsolve(l, x, k=ncol(l), upper.tri=FALSE, transpose=FALSE)
这两个函数都是符合操作的函数，大致可以分为三个步骤：
①通过将系数矩阵的上三角或者下三角变为0的到新的系数矩阵,这通过upper.tri参数来实现，若upper.tri=TRUR,上三角不为0。
②通过将对步骤1中得到的新系数矩阵进行转置得到新的系数矩阵，这通过transpose参数实现，若transpose=FALSE，则步骤1中得到的系数矩阵将被转置。
③根据步骤2得到的系数矩阵来解方程组。
X1+4X2+7X3=1
2X1+5X2+8X3=2
3X1+6X2+9X3=3
方程组的系数矩阵为：
> A
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9
> b
[1] 1 2 3
> backsolve(A,b,upper.tri=T,transpose=F)
[1] -0.8000000 -0.1333333 0.3333333
过程分解：
①upper.tri=T，说明系数矩阵的上三角不为0。
> B=A
> B[lower.tri(B)]=0
> B
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 0 5 8
[3,] 0 0 9
②transpose=F说明系数矩阵未被转置。
③解方程：
> solve(B,b)
[1] -0.8000000 -0.1333333 0.3333333
5 其它
5.1矩阵的向量化
将矩阵向量化有时候是必要的。矩阵的向量化可以通过as.vector()来实现：
> A
[,1] [,2] [,3] [,4]
[1,] 1 4 7 10
[2,] 2 5 8 11
[3,] 3 6 9 12
将矩阵元素向量化：
> as.vector(A)
[1] 1 2 3 4 5 6 7 8 9 10 11 12
将矩阵的方阵部分元素向量化：
> as.vector(A[1:min(dim(A)),1:min(dim(A))])
[1] 1 2 3 4 5 6 7 8 9
5.2矩阵的合并
5.2.1矩阵的列合并
矩阵的列合并可以通过cbind()来实现。
> A
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9
> B=1:3
> cbind(A,B)
B
[1,] 1 4 7 1
[2,] 2 5 8 2
[3,] 3 6 9 3
5.2.2矩阵的行合并
矩阵的行合并可以通过rbind()来实现。
> A
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9
> B=1:3
> rbind(A,B)
[,1] [,2] [,3]
1 4 7
2 5 8
3 6 9
B 1 2 3
5.3 时序矩阵的滞后
在时间序列中经常会用到一个序列的滞后序列，R中的包fMultivar中的函数tslag()提供了这个功能。
> library(fMultivar)
Loading required package: sn
Loading required package: mnormt
Package ‘sn’, 0.4-16 (2010-08-30). Type ‘help(SN)’ for summary information
Loading required package: timeDate
Loading required package: timeSeries
Loading required package: fBasics
Loading required package: MASS
Attaching package: ‘fBasics’
The following object(s) are masked from ‘package:base’:
norm
> args(tslag)
function (x, k = 1, trim = FALSE)
NULL
其中：x为一个向量，k指定滞后阶数，可以是一个自然数列，若trim为假，则返回序列与原序列长度相同，但含有NA值；若trim项为真，则返回序列中不含有NA值，例如：
> x=1:9
> x
[1] 1 2 3 4 5 6 7 8 9
> 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
> 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
=============================
主要包括以下内容：
创建矩阵向量；矩阵加减，乘积；矩阵的逆；行列式的值；特征值与特征向量；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为列数，注意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
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
> 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,] 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
> solve (a) %*%a
[,1] [,2] [,3] [,4]
[1,] 1.000000e+00 2.748453e-17 -2.787755e-17 -8.023096e-17
[2,] 1.626303e-19 1.000000e+00 -4.960225e-18 6.977925e-16
[3,] 2.135878e-17 -4.629543e-17 1.000000e+00 6.201636e-17
[4,] 1.866183e-17 1.563962e-17 1.183813e-17 1.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,] 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
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)
> 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
对于函数forwardsolve()而言，
例如：
> A
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 2 5 8
[3,] 3 6 9
> B
[,1] [,2] [,3]
[1,] 1 0 0
[2,] 2 5 0
[3,] 3 6 9
> C
[,1] [,2] [,3]
[1,] 1 4 7
[2,] 0 5 8
[3,] 0 0 9
> x
[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
这两个函数同样可以用于取一个矩阵的上下三角矩阵，例如：
> 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
21 行列式的值
在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
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)
> x
[,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本文章为转载内容，我们尊重原作者对文章享有的著作权。如有内容错误或侵权问题，欢迎原作者联系我们进行内容更正或删除文章。
