R語言－強大的矩陣運算

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維矩陣，公式爲：

              a^­₁₁B … a_1nB

A_m×n×B_h×k=     …     …

               a_m1B … a_mnB   _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爲方程組的右端。例如：

已知方程組：

2x₁+2x₃=1

2x₁+x₂₊2x­₃=2

2x₁+x­₂=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，x₂=1，x₃=-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獲得的係數矩陣來解方程組。

X₁+4X₂+7X₃=1

2X₁+5X₂+8X₃=2

3X₁+6X₂+9X₃=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