# 【3.2】矩阵操作（终结版）

## 一  矩阵基本操作

### 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


> X=seq(1,4,length=4)
> X
[1] 1 2 3 4


### 1.2 创建矩阵

R中创建矩阵的方法也有很多。大致分为直接创建和由其它格式转换两种方法。

### 1.2.1直接创建矩阵

> args(matrix)
function (data = NA, nrow = 1, ncol = 1, byrow = FALSE, dimnames = NULL)
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


> 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 查看和改变矩阵的维数

> 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(X)
[1] 4

> dim(X)[2]
[1] 3

> ncol(X)
[1] 3

> 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 矩阵行列的名称

> 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

> 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

> 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

> 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.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 矩阵的行和与列和及行平均值和列均值

> 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 矩阵的数乘

> 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



### 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

> 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


### 2.6.2 矩阵的Kronecker积

n×m矩阵A和h×k矩阵B的Kronecker积是一个nh×mk维矩阵，公式为：

 a­11B … a1nB
Am×n×Bh×k= … …
am1B … amnB mh×nk

> 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 矩阵的伴随矩阵

>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
[,1] [,2] [,3]
[1,] -3 6 -3
[2,] 6 -12 6
[3,] -3 6 -3


### 2.8.1 矩阵的逆

> 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

> 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）

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



### 2.8.3 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矩阵的特征值和特征向量

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

> 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

> 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.1 Choleskey分解

> 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

> t(chol(A))%*%chol(A)
[,1] [,2] [,3]
[1,] 2 1 1
[2,] 1 2 1
[3,] 1 1 2

> crossprod(chol(A),chol(A))
[,1] [,2] [,3]
[1,] 2 1 1
[2,] 1 2 1
[3,] 1 1 2

> prod(diag(chol(A))^2)
[1] 4
> det(A)
[1] 4

> 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.1 普通方程组

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



### 4.2 特殊方程组

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


1. 通过将系数矩阵的上三角或者下三角变为0的到新的系数矩阵,这通过upper.tri参数来实现，若upper.tri=TRUR,上三角不为0。
2. 通过将对步骤1中得到的新系数矩阵进行转置得到新的系数矩阵，这通过transpose参数实现，若transpose=FALSE，则步骤1中得到的系数矩阵将被转置。
3. 根据步骤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.1矩阵的向量化

> 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.1矩阵的列合并

> 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矩阵的行合并

> 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 时序矩阵的滞后

> library(fMultivar)
Package 'sn', 0.4-16 (2010-08-30). Type 'help(SN)' for summary information
Attaching package: 'fBasics'
The following object(s) are masked from 'package:base':
norm
> args(tslag)
function (x, k = 1, trim = FALSE)
NULL

> 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


http://blog.renren.com/share/159280283/14363739061