R语言矩阵运算 – Neohope的网络笔记

1向量
1.1定义向量

##定义向量，下面两种方式结果相同
##c(..., recursive = FALSE) 
x=1:12
x=c(1,2,3,4,5,6,7,b,9.10)

1.2向量与数字四则运算

##+-×/
x=1:12
y=10
z=x*y

1.3向量的内积

##%*% crossprod
x=1:5
y=1:5
z=x%*%y
z1=crossprod(x)
z2=crossprod(x,y)

1.4向量的外积

##x%o%y tcrossprod
x=1:5
y=1:5
z=x%o%y
z1=tcrossprod(x)
z2=tcrossprod(x,y)
z3=outer(x,y)

1.5矩阵转向量

vec<-function (x){
          t(t(as.vector(x)))
}
vech<-function (x){
          t(x&#91;lower.tri(x,diag=T)&#93;)
}
x=matrix(1:16,nrow=4,ncol=4,byrow=TRUE)
vec(x)
vech(x)
&#91;/code&#93;

1.6生成滞后序列
&#91;code lang="R"&#93;
library(fMultivar)
x=1:20
tslag(x,1:4,trim=F)
tslag(x,1:4,trim=T)
&#91;/code&#93;

<strong>2数组</strong>
2.1定义数组
[code lang="R"]
##向量转数组，按列转
##dim(x)
dim(x)=c(3,4);
##或者直接定义多维数组
x=array(1:12,dim=c(3,4))

3矩阵
3.1定义矩阵

##matrix(data=NA，nrow=1，ncol=1，byrow=FALSE，dimnames=NULL)
##data项为必要的矩阵元素，nrow为行数，ncol为列数，nrow与ncol的乘积应为矩阵元素个数，byrow项控制排列元素时是否按行进行，dimnames给定行和列的名称
x=matrix(1:16,nrow=4,ncol=4,byrow=TRUE)
##行数
nrow(x)
##列数
ncol(x)
##上三角阵
C=x
C[lower.tri(A)]=0
##下三角阵
C=x
C[upper.tri(A)]=0
##元素的行与列
row(x)
col(x)
D=x
D[row(x)<col(x)&#93;=0
D=x
D&#91;row(x)>col(x)]=0

3.2定义对角矩阵

##diag
A=c(1:5)
##对角阵
B=diag(A)
##获取对角元素
C=diag(B)
##生成单位对角阵
D=diag(5)

3.3矩阵与数字四则运算

##+-*/
A=matrix(1:9,nrow=3,ncol=3,byrow=TRUE)
B=10
C=A+B
D=A-B
E=A*B
F=A/B
##行和
rowSums(A)
apply(A,1,sum)
##行平均值
rowMeans(A)
apply(A,1,mean)
##行方差
apply(A,1,var)
##列和
colSums(A)
apply(A,2,sum)
##列平均值
colMeans(A)
apply(A,2,mean)
##列方差
apply(A,2,var)

3.4矩阵与矩阵加减法

##+-
A=matrix(1:9,nrow=3,ncol=3,byrow=TRUE)
B=matrix(10:18,nrow=3,ncol=3,byrow=TRUE)
C=A+B
D=A-B

3.5矩阵与矩阵乘法

##%*%
A=matrix(1:9,nrow=3,ncol=3,byrow=TRUE)
B=matrix(10:18,nrow=3,ncol=3,byrow=TRUE)
C=A%*%B
##t(A)%*%B;
C1=crossprod(A,B)
##A%*%t(B)
C2=tcrossprod(A,B)
##Hadamard积
D=A*A

3.6矩阵转置

##t(x)
A=matrix(1:6,nrow=2);
B=t(A)

3.7矩阵行列式

##det(x, ...)
A=matrix(1:9,nrow=3);
B=det(A)

3.8矩阵求逆

##solve
A=matrix(c(1,2,3,2,2,1,3,4,3),nrow=3,byrow=TRUE);
B=det(A)
C=solve(A)

##需要strucchange包,A'A的逆
library(strucchange)
C1=solveCrossprod(A,method="qr")
C2=solveCrossprod(A,method="chol")
C3=solveCrossprod(A,method="solve")
C4=solve(crossprod(A,A))

3.9求解线性方程

##solve
A=matrix(c(1,2,3,2,2,1,3,4,3),nrow=3,byrow=TRUE);
B=c(1,2,8)
C=solve(A,B)

##对于系数为三角阵的情况，可以用 backsolve及fowardsolve函数求解
##这一段有些错误，时间不够了，下次修复
A=matrix(1:9,3,3)
V=c(1,2,3)
B=A
B[upper.tri(B)]=0
backsolve(A,V,upper.tri=F,transpose=F)
forwardsolve(A,V,upper.tri=F,transpose=F)
##why not the same?
##something is wrong
solve(B,V)

3.10特征值及特征向量

##eigen
A=matrix(c(1,-2,2,-2,-2,4,2,4,-2),nrow=3,byrow=TRUE);
B=eigen(A)
##结果应该为对角阵
C=solve(B$vectors)%*%A%*%B$vectors

4.矩阵分解
4.1奇异值分解

##A=U%*%D%*%t(V)，其中U, V是正交阵，D为对角阵，也就是矩阵A的奇异值
##svd
A=matrix(c(1,0,1,0,0,1),nrow=3,byrow=TRUE);
B=svd(A)
C=B$u%*%diag(B$d)%*%t(B$v)

4.2QR分解

##设A为m*n矩阵，如果存在m*m酉矩阵Q（即Q(H)Q=QQ(H)=I）和m*n阶梯形矩阵R，使得A=QR，那么此分解称为QR分解
##qr
A=matrix(c(1,1,1,2,-1,-1,2,-4,5),nrow=3,byrow=TRUE);
B=qr(A)
Q=qr.Q(B)
R=qr.R(B)
C=Q%*%R

4.3Schur分解

##A=USU(H)，其中U是标准的正交矩阵（即满足UU(H)=I），S为上三角矩阵，并且对角线上的元素为A的特征值。
##Schur
library(Matrix);
A=matrix(c(1,1,1,2,-1,-1,2,-4,5),nrow=3,byrow=TRUE);
B=Schur(A, vectors=TRUE);
C=B$Q%*%B$T%*%t(B$Q)

4.4Cholesky分解

##对任意的正定矩阵A，存在上三角矩阵R，使A=t(R)%*%R，则称为A的Cholesky分解
##chol
A=matrix(c(4,-1,1,-1,4.25,2.75,1,2.75,3.5),nrow=3,byrow=TRUE);
B=chol(A);
C=t(B)%*%B;
##求行列式
D=prod(diag(chol(A))^2)
##求逆
E=chol2inv(chol(A))

Leave a Reply Cancel reply