R语言矩阵运算

爱被打了一巴掌 2022-08-12 04:29 827阅读 0赞

# R语言矩阵运算    主要包括以下内容：   创建矩阵向量；矩阵加减，乘积；矩阵的逆；行列式的值；特征值与特征向量；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