這是我見過最全的R矩陣運(yùn)算的例子。
原文：
https://blog.csdn.net/kkkkkiko/article/details/83118783
PS：這位作者轉(zhuǎn)的其他人的，但是原文已經(jīng)不見了。以防這個(gè)也消失了。。。。本人再轉(zhuǎn)一次，以防丟失這么好的例子。

1 矩陣基本操作

1.1創(chuàng)建向量

R里面有多種方法來創(chuàng)建向量（Vector），最簡(jiǎn)單的是用函數(shù)c()。例如：

>X=c(1,2,3,4)

>X

[1] 1 2 3 4

當(dāng)然，還有別的方法。例如：

>X=1:4

>X

[1] 1 2 3 4

還有seq()函數(shù)。例如：

> X=seq(1,4,length=4)

> X

[1] 1 2 3 4

注意一點(diǎn)，R中的向量默認(rèn)為列向量，如果要得到行向量需要對(duì)其進(jìn)行轉(zhuǎn)置。

1.2創(chuàng)建矩陣

R中創(chuàng)建矩陣的方法也有很多。大致分為直接創(chuàng)建和由其它格式轉(zhuǎn)換兩種方法。

1.2.1直接創(chuàng)建矩陣

最簡(jiǎn)單的直接創(chuàng)建矩陣的方法是用matrix()函數(shù)，matrix()函數(shù)的使用方法如下：

> args(matrix)

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

NULL

其中，data參數(shù)輸入的為矩陣的元素，不能為空；nrow參數(shù)輸入的是矩陣的行數(shù)，默認(rèn)為1；ncol參數(shù)輸入的是矩陣的列數(shù)，默認(rèn)為1；byrow參數(shù)控制矩陣元素的排列方式，TRUE表示按行排列，FALSE表示按列排列，默認(rèn)為FALSE；dimnames參數(shù)輸入矩陣的行名和列名，可以不輸入，系統(tǒng)默認(rèn)為NULL。例如：

> matrix(1:6,nrow=2,ncol=3,byrow=FALSE)

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

[1,] 1 3 5

[2,] 2 4 6

改變矩陣的行數(shù)和列數(shù)：

> matrix(1:6,nrow=3,ncol=2,byrow=FALSE)

** [,1] [,2]**

[1,] 1 4

[2,] 2 5

[3,] 3 6

改變byrow參數(shù)：

> matrix(1:6,nrow=3,ncol=2,byrow=T)

** [,1] [,2]**

[1,] 1 2

[2,] 3 4

[3,] 5 6

設(shè)定矩陣的行名和列名：

> 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 由其它格式轉(zhuǎn)換

也可以由其它格式的數(shù)據(jù)轉(zhuǎn)換為矩陣，此時(shí)需要用到函數(shù)as.matrix()。

1.3 查看和改變矩陣的維數(shù)

矩陣有兩個(gè)維數(shù)，即行維數(shù)和列維數(shù)。在R中查看矩陣的行維數(shù)和列維數(shù)可以用函數(shù)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

只返回行維數(shù)：

> dim(X)[1]

[1] 4

也可以用函數(shù)nrow()

> nrow(X)

[1] 4

只返回列維數(shù)：

> dim(X)[2]

[1] 3

也可以用函數(shù)ncol():

> ncol(X)

[1] 3

同時(shí)，函數(shù)dim()也可以改變矩陣的維數(shù)。例如：

> 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矩陣行列的名稱

查看矩陣的行名和列名分別用函數(shù)rownames()和函數(shù)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”

同時(shí)也可以改變矩陣的行名和列名，比如：

> 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 矩陣元素的查看及其重新賦值

查看矩陣的某個(gè)特定元素，只需要知道該元素的行坐標(biāo)和列坐標(biāo)即可，例如：

> 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

也可以重新對(duì)矩陣的元素進(jìn)行賦值，將矩陣第二行、第三列的元素替換為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中有一個(gè)diag()函數(shù)可以返回矩陣的全部對(duì)角元素：

> 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

當(dāng)然也可以對(duì)對(duì)角元素進(jìn)行重新賦值：

> diag(X)=c(0,0,1)

> X

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

[1,] 0 4 7

[2,] 2 0 8

[3,] 3 6 1

當(dāng)操作對(duì)象不是對(duì)稱矩陣時(shí)，diag()也可以進(jìn)行操作。

> 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()函數(shù)還能用來生成對(duì)角矩陣：

> diag(c(1,2,3))

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

[1,] 1 0 0

[2,] 0 2 0

[3,] 0 0 3

也可以生成單位對(duì)角矩陣：

> 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中查看矩陣的上三角和下三角部分很簡(jiǎn)單?？梢酝ㄟ^lower.tri()和upper.tir()來實(shí)現(xiàn)：

> 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 矩陣計(jì)算

2.1矩陣轉(zhuǎn)置

R中矩陣的轉(zhuǎn)置可以用t()函數(shù)完成，例如：

> 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中很容易計(jì)算一個(gè)矩陣的各行和和各列和以及各行的平均值和各列的平均值。例如：

> 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中的函數(shù)det()將計(jì)算方陣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矩陣相加減

矩陣元素的相加減是指維數(shù)相同的矩陣，處于同行和同列的位置的元素進(jìn)行加減。實(shí)現(xiàn)這個(gè)功能用“＋”，“－”即可。例如：

> 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矩陣的數(shù)乘

矩陣的數(shù)乘是指一個(gè)常數(shù)與一個(gè)矩陣相乘。設(shè)A為m****×****n矩陣，c****≠****0，在R中求cA的值，可以用符號(hào)“*”。例如：

> 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

> cA*

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

[1,] 2 8 14 20

[2,] 4 10 16 22

[3,] 6 12 18 24

結(jié)果矩陣與原矩陣的所有相應(yīng)元素都差一個(gè)常數(shù)c。

2.6矩陣相乘

2.6.1矩陣的乘法

A****為****m****×****n矩陣，B為n******×******k矩陣，在R中求AB，可以符號(hào)“%%****”*。例如：

> 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無意義，因其不符合矩陣的相乘規(guī)則。

若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

也可以用函數(shù)crossprod()計(jì)算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積是一個(gè)nh****×mk維矩陣，公式為：

** a-11B ****… a1nB**

Am****×n****×****Bh****×k= ****… ****…

** am1B ****… amnB mh****×nk**

在R中Kronecker積可以用函數(shù)kronecher()來計(jì)算。例如：

> 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包中的函數(shù)make.adjoint()函數(shù)。例如：

>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可以用函數(shù)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

驗(yàn)證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函數(shù)可以更好的得到結(jié)果：

> round(A%%solve(A))*

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

[1,] 1 0 0

[2,] 0 1 0

[3,] 0 0 1

solve()函數(shù)也可以用來求解方程組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()函數(shù)可以計(jì)算矩陣的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

驗(yàn)證性質(zhì)①：

> 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

驗(yàn)證性質(zhì)②：

> 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

驗(yàn)證性質(zhì)③：

> 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

驗(yàn)證性質(zhì)④：

> 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

也可以不必如此麻煩來驗(yàn)證性質(zhì)③和④，因?yàn)棰酆廷苤皇潜砻?strong>AA+和A+A是對(duì)稱矩陣。

2.8.3 X’X的逆

很多時(shí)候，我們需要計(jì)算形如X’X的逆。這很容易實(shí)現(xiàn)，例如：

> 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包中的函數(shù)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矩陣的特征值和特征向量

可以通過對(duì)矩陣A進(jìn)行譜分解來得到矩陣的特征值和特征向量。矩陣A的譜分解如下：A=U****Λ****U’，其中U的列為A的特征值所對(duì)應(yīng)的特征向量，在R中可以用eigen()函數(shù)得到U和Λ。例如：

> args(eigen)

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

NULL

其中，x參數(shù)輸入矩陣；symmetric參數(shù)判斷矩陣是否為對(duì)稱矩陣，如果參數(shù)為空，系統(tǒng)將自動(dòng)檢測(cè)矩陣的對(duì)稱性。例如：

> 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 矩陣高級(jí)操作

3.1 Choleskey分解

對(duì)于正定矩陣A，可以對(duì)其進(jìn)行Choleskey分解，A=P’P，其中P為上三角矩陣，在R中可以用函數(shù)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

驗(yàn)證A=P’P****：

> t(chol(A))%%chol(A)*

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

[1,] 2 1 1

[2,] 1 2 1

[3,] 1 1 2

也可以用crossprod()函數(shù)：

> crossprod(chol(A),chol(A))

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

[1,] 2 1 1

[2,] 1 2 1

[3,] 1 1 2

可以用Choleskey分解來計(jì)算矩陣的行列式：

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

[1] 4

> det(A)

[1] 4

也可以用Choleskey分解來計(jì)算矩陣的逆，這時(shí)候可以用到函數(shù)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中可以用函數(shù)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矩陣可以進(jìn)行QR分解:A=QR，其中Q’Q=I，在R中可以用函數(shù)qr()來完成這個(gè)過程，例如：

> 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項(xiàng)包含了Q矩陣和R矩陣的信息。要想得到Q矩陣和R矩陣，可以用qr.Q()函數(shù)和qr.R()函數(shù)：

> 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普通方程組

解普通方程組可以用函數(shù)solve()，solve()的基本用法是solve(A,b)，其中，A為方程組的系數(shù)矩陣，b為方程組的右端。例如：

已知方程組：

2x1+2x3=1

2x1+x2+2x-3=2

2x1+x-2=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

即x-1=1，x2=1，x3=-0.5。

4.2 特殊方程組

對(duì)于系數(shù)矩陣是上三角矩陣和下三角矩陣的方程組。R中提供了backsolve()和fowardsolve()來解決這個(gè)問題。

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

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

這兩個(gè)函數(shù)都是符合操作的函數(shù)，大致可以分為三個(gè)步驟：

①通過將系數(shù)矩陣的上三角或者下三角變?yōu)?的到新的系數(shù)矩陣,這通過upper.tri參數(shù)來實(shí)現(xiàn)，若upper.tri=TRUR,上三角不為0。

②通過將對(duì)步驟1中得到的新系數(shù)矩陣進(jìn)行轉(zhuǎn)置得到新的系數(shù)矩陣，這通過transpose參數(shù)實(shí)現(xiàn)，若transpose=FALSE，則步驟1中得到的系數(shù)矩陣將被轉(zhuǎn)置。

③根據(jù)步驟2得到的系數(shù)矩陣來解方程組。

X1+4X2+7X3=1

2X1+5X2+8X3=2

3X1+6X2+9X3=3

方程組的系數(shù)矩陣為：

> 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，說明系數(shù)矩陣的上三角不為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說明系數(shù)矩陣未被轉(zhuǎn)置。

③解方程：

> solve(B,b)

[1] -0.8000000 -0.1333333 0.3333333

5 其它

5.1矩陣的向量化

將矩陣向量化有時(shí)候是必要的。矩陣的向量化可以通過as.vector()來實(shí)現(xiàn)：

> 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()來實(shí)現(xiàn)。

> 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()來實(shí)現(xiàn)。

> 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 時(shí)序矩陣的滯后

在時(shí)間序列中經(jīng)常會(huì)用到一個(gè)序列的滯后序列，R中的包fMultivar中的函數(shù)tslag()提供了這個(gè)功能。

> 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為一個(gè)向量，k指定滯后階數(shù)，可以是一個(gè)自然數(shù)列，若trim為假，則返回序列與原序列長(zhǎng)度相同，但含有NA值；若trim項(xiàng)為真，則返回序列中不含有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