R语言matrix运算

### 1. 矩阵创建
# M <- matrix(1:12,3,4)
# args(matrix)
M <- matrix(sample(0:1000,1000,replace=TRUE),nrow=100,ncol=10,
            dimnames=list(paste0('g',1:100),paste0('s',1:10)))
nrow(M)
ncol(M)
rownames(M)
colnames(M)

A=1:15
class(A) # [1] "integer"
dim(A)=c(3,5)
class(A) # [1] "matrix" "array" 

### 2. 矩阵转置
M <- matrix(1:15,nrow=3,ncol=5)
t(M) # 转置，行列互换

## 向量转置后变矩阵
x <- 1:15
X <- t(x) # 只有一列的矩阵
class(X)  # [1] "matrix" "array"

class(t(X)) # 还是矩阵
### 3.矩阵求逆
M <- matrix(rnorm(9),nrow=3,ncol=3) 

# solve(a, b)运算结果是解线性方程组ax = b，
# 若b缺省，则系统默认为单位矩阵，因此可用其进行矩阵求逆
solve(M) # 要求是方阵

solve(M) %*% M   # %*%  矩阵乘法运算符

### 4. 矩阵元素加减乘除
M1 <- matrix(rnorm(9,10,1),nrow=3,ncol=3)
M2 <- matrix(rnorm(9,10,1),nrow=3,ncol=3)
M1 + M2; M1 - M2; M1 / M2;
round(M1 * M2, 2) #  对应元素相乘
M1 + 6; M1 - 6; M1 * 6; M1 / 6;

### 5. 矩阵相乘
A <- matrix(1:15,nrow=3,ncol=5)
B <- matrix(1:15,nrow=5,ncol=3)

A %*% B # A矩阵列数等于B矩阵的行数

### 6. 内积和外积
A1 <- matrix(rnorm(15),nrow=3,ncol=5)
A2 <- matrix(rnorm(15),nrow=3,ncol=5)
#dim(A1) # [1] 3 5
#dim(A2) # [1] 3 5
A3 <- crossprod(A1,A2) # 内积
#dim(A3) # [1] 5 5
A4 <- t(A1) %*% A2
identical(A3,A4)

tcrossprod(A1,A2) # 外积
identical(tcrossprod(A1,A2),A1 %*% t(A2))

(z <- crossprod(1:4))    # = sum(1 + 2^2 + 3^2 + 4^2)
drop(z)                  # 删除只有一个级别的数组的维度。scalar
x <- 1:4; names(x) <- letters[1:4]; x
tcrossprod(as.matrix(x)) # is
identical(tcrossprod(as.matrix(x)),
          crossprod(t(x)))
tcrossprod(x)            # no dimnames

m <- matrix(1:6, 2,3) ; v <- 1:3; v2 <- 2:1
stopifnot(identical(tcrossprod(v, m), v %*% t(m)),
          identical(tcrossprod(v, m), crossprod(v, t(m))),
          identical(crossprod(m, v2), t(m) %*% v2))

### 7.行列统计
A <- matrix(rnorm(15),nrow=3,ncol=5)
rowSums(A)
rowMeans(A)
colSums(A)
colMeans(A)

normalize <- function(x){
    return ((x-min(x))/(max(x)-min(x)))
}

apply(A,2,normalize) # 按列进行处理


### 8.对角阵
# 提取或替换矩阵的对角线，或构造对角线矩阵。
diag(1,nrow=3,ncol=3)
diag(c(5,9,12,6),nrow=4,ncol=4)
A <- matrix(1:9,nrow=3,ncol=3)
diag(A) # 提取A的对角线元素
diag(diag(A))


### 9. 子矩阵
A <- matrix(1:15,nrow=3,ncol=5)
A[1:2,1:3]
#取矩阵的上、下三角部分
lower.tri(A) # 返回TRUE FALSE矩阵
lower.tri(A,diag=T) # 加上对角线上的元素
A[lower.tri(A)]=0

upper.tri(A)
upper.tri(A,diag=T)
A[upper.tri(A,diag=T)]=0

### 10.矩阵的特征值与特征向量
# 矩阵A的谱分解为A=UΛU‘,其中Λ是由A的特征值组成的对角矩阵，
# U的列为A的特征值对应的特征向量
## 要求是方阵
M1 <- matrix(1:25,nrow=5,ncol=5)
eigen1=eigen(M1,symmetric=T)
eigen2=eigen(M1)

### 11.矩阵广义逆
library("MASS")
M1 <- matrix(1:15,nrow=3,ncol=5) 
dim(M1) #[1] 3 5
M2 <- ginv(M1)
dim(M2) #[1] 5 3

identical(M1%*%ginv(M1)%*%M1,M1) #[1] FALSE
# M1%*%ginv(M1)%*%M1 近似等于M1

identical(ginv(M1)%*%M1%*%ginv(M1),ginv(M1)) #[1] FALSE
# ginv(M1)%*%M1%*%ginv(M1) 近似等于ginv(M1)

identical(t(M1%*%ginv(M1)),M1%*%ginv(M1)) #[1] FALSE
# t(M1%*%ginv(M1))近似等于M1%*%ginv(M1)

identical(t(ginv(M1)%*%M1),ginv(M1)%*%M1) #[1] FALSE
# t(ginv(M1)%*%M1)近似等于 ginv(M1)%*%M1

### 11.矩阵奇异值分解
M <- matrix(1:15,nrow=3,ncol=5) 
svd(M)

### 12.矩阵QR分解
M <- matrix(1:15,nrow=3,ncol=5) 
qr(M)
qr.R(qr(M)) # 矩阵R
qr.Q(qr(M)) # 矩阵Q

### 13.矩阵行列式
# 必需是方阵
M <- matrix(rnorm(25),nrow=5,ncol=5) 
det(M)

(x <- matrix(1:4, ncol = 2))
unlist(determinant(x))
det(x)

det(print(cbind(1, 1:3, c(2,0,1))))