1.3.3 矩阵

1.3.3.1 矩阵结构

矩阵（matrix）是 n 行和 p 列数值的长方形排列，是二维数组（two－dimensional array）。矩阵 A （ n × p ）的形式为：

矩阵中每个元素（如 a ₁₁ 和 a _{n 2} ）的数据类型相同（如数值型、字符型或逻辑型）。R内存矩阵的方式体现带有维度属性（dimension attribute）的数值向量，内存数据框的方式是数据列。矩阵的维度指矩阵中的行数和列数。

假如有 n ＝20名英语学习者在三个因变量（dv1、dv2和dv3）上被测量，数据存储较为方便的方式是利用函数cbind创建一个由 n ＝20排和 p ＝3列构成的矩阵。例如：

＞dv1＜－c（60，65，70，70，75，78，80，63，90，67，66，77，84，80，90，75，88，63，82，80）

＞dv2＜－c（66，67，77，72，60，67，78，65，87，68，60，70，73，70，85，62，80，60，70，66）

＞dv3＜－c（61，62，66，65，70，69，81，60，85，62，60，70，75，73，80，70，77，60，71，69）

＞z＜－cbind（dv1，dv2，dv3）

＞z［1：5，1：3］

dv1 dv2 dv3

［1，］ 60 66 61

［2，］ 65 67 62

［3，］ 70 77 66

［4，］ 70 72 65

［5，］ 75 60 70

由于数据集较大，这里只显示3列前5行的数据。另外一种创建矩阵的方式是利用函数matrix。例如：

＞zz＜－matrix（c（dv1，dv2，dv3），ncol＝3）

＞zz［1：5，1：3］

［，1］［，2］［，3］

［1，］ 60 66 61

［2，］ 65 67 62

［3，］ 70 77 66

［4，］ 70 72 65

［5，］ 75 60 70

矩阵z和zz相同，只是矩阵z显示变量名，矩阵zz显示变量所在的列。

要确定矩阵包括的行数和列数，可以利用函数nrow和ncol，如nrow（zz）。对矩阵的行和列可以进行算术运算。譬如，利用mean（z［，1］）计算矩阵z第一列的平均数，得到75.15。如果要计算所有行和列的平均数或其他统计量，一个简便的方法是使用函数apply。函数apply既适用于矩阵又适用于数据框中行和列向量的计算。譬如，执行命令apply（z，2，mean）得到每列向量的平均数：

＞apply（z，2，mean）＃for a matrix，2 indicates columns

dv1 dv2 dv3

75.15 70.15 69.30

执行命令apply（z，1，mean）得到每行向量的平均数：

＞apply（z，1，mean）＃for a matrix，1 indicates rows

［1］62.33333 64.66667 71.00000 69.00000 68.33333 71.33333 79.66667 62.66667

［9］87.33333 65.66667 62.00000 72.33333 77.33333 74.33333 85.00000 69.00000

［17］81.66667 61.00000 74.33333 71.66667

1.3.3.2 特殊矩阵

在统计分析中会用到一些特殊的矩阵。本节简要介绍方阵、转置矩阵、对角矩阵、单位阵和正交矩阵。

方阵（square matrix）。 行数和列数相同（ r ＝ c ）的矩阵是方阵一。种常见的方阵是相关矩阵，记作 R ＝（ r _ij ）。例如：

在这个矩阵中，第一个变量与第二个变量的相关系数是0.2，第一个变量与第三个变量的相关系数是0.5，第二个变量与第三个变量的相关系数是0.3。对角线上的值为1，因为变量与自身完全相关。矩阵是对称的。

转置矩阵 （ transpose of a matrix ）。一个矩阵的行变成相应的列得到的矩阵为转置矩阵。矩阵 X 的转置矩阵记作 X′ 或 X ^T 。例如：

相关矩阵 R 的转置矩阵 R ′与矩阵 R 相同，因为相关矩阵是对称矩阵。

对角矩阵 （ diagonal matrix ）。对角线之外所有元素的值为0的方阵是对角矩阵。例如：

对角矩阵R函数是diag。上例对角矩阵的R命令为diag（c（3，5，8））。

单位阵 （ identity matrix ， I ）。如果对角矩阵对角线上元素的值均为1，则为单位阵，记作 I 。例如：

单位阵的R函数也是diag。上例单位阵可用R命令diag（3）或diag（rep（1，3））得到。

正交矩阵 （ orthogonal matrix ）。正交矩阵同对角矩阵一样也是一种特殊形式的方阵。如果一个 p × p 方阵 Q 满足条件式 QQ ′＝ Q ′ Q ＝ I _p ，方阵 Q 则为正交矩阵。例如：

正交矩阵的一个特点是，如果矩阵 Q 有一个 i 排 表明 ，其中 i ≠ j 。譬如，前面的矩阵 A 有以下属性：

1.3.3.3 矩阵计算

矩阵可以进行加、减、乘、除计算。本小节利用数值例子介绍矩阵计算方法。

加减法 。如果矩阵与一个标量进行加、减、乘、除运算，则按矩阵元素计算的方法，即对矩阵的每个元素进行加、减、乘、除运算。譬如， 。对于包括标量的矩阵计算，R使用加、减、乘、除的计算符依次为：＋、－、＊、／。例如：

两个矩阵加、减时，矩阵的维度（即排数和列数）要相同。计算时，矩阵对应元素相加或相减。例如：

点乘 （ dot product ）。矩阵点乘时，一个矩阵的列数必须等于另一个矩阵的排数。计算时，第一个矩阵第一排上的元素乘以另一个矩阵第一列上的对应元素，然后将乘积相加，得到第一排和第一列位置上的值；第一个矩阵第一排上的元素乘以另一个矩阵第二列上的对应元素，然后将乘积相加，得到第一排和第二列位置上的值，以此类推。例如：

本例中，第一个矩阵为2×3矩阵，第二个矩阵是3×3方阵，第一个矩阵的列数等于第二个矩阵的排数，得到的点乘结果为2×3矩阵。矩阵点乘必须满足第一个矩阵的列数等于第二个矩阵的排数这一条件。能够相乘的两个矩阵是相容的（conformable）。

R使用的点乘符号是％＊％。注意，两个矩阵位置不同，点乘结果可能不同。例如，矩阵 A 与 C 点乘的结果为：

为了增加对矩阵点乘的了解，我们使用前面的for循环编写以下R函数：

Mat.mult＝function（A，B）｛

m＝nrow（A）

n＝ncol（A）

p＝ncol（B）

C＝matrix（0，nrow＝m，ncol＝p）

for（i in 1：m）｛

for（j in 1：p）｛

sumvalue＝0

for（k in 1：n）｛

sumvalue＝sumvalue＋A［i，k］＊B［k，j］

｝

C［i，j］＝sumvalue

｝

｝

return（C）

｝

利用编写的函数Mat.mult计算矩阵 A 与 C 点乘的R命令和计算结果为：

以上结果表明，利用函数Mat.mult的计算结果与利用点乘符号％＊％的计算结果相同。

点乘是非常有用的统计量计算方法。如果随机变量 X 是包括各个观测值 X ₁ ，…， 的列向量， J 是排向量（1，…，1）（ n 个 1），那么样本平均数就为： ，方差（variance， ）就为：

克罗内克积 （ Kronecker product ）。如果 A 是 m × n 矩阵， B 是 p × q 矩阵，克罗内克积 是 mp × nq 矩阵，形式为：

例如：

已知 D 矩阵为D＜－matrix（c（5，6，4，7，3，2，2，1，8），ncol＝3），则 是9×9矩阵。R计算克罗内克积的函数是kronecker，如kronecker（C，D）。

方阵的迹（trace）。 一个方阵对角元素之和是方阵的迹，常记作tr。例如，若 D ＝ ，则tr（D）＝5＋3＋8＝16。R没有直接的函数计算方阵的迹，但是可以很容易地用sum（diag（））命令得到。我们也可以调用R数据包psych中的函数tr（m），其中m是方阵。

行列式 （ determinant ）。行列式是与任一方阵 A 相联系的一个值，记作｜ A ｜。行列式是实数（real number），仅存在于方阵中。对于一个2×2矩阵 A ，如 A ＝ ，行列式为｜ A ｜＝ ad － bc ，即主对角线元素乘积与主对角线外元素乘积之差。行列式的一个重要特点是，如果方阵中至少有两个向量存在线性相依（linearly dependent），则行列式为0。通常，随着方阵中各个向量独立性的增强，行列式越偏离0值。

如果方阵是高阶矩阵，行列式计算比较复杂。这里举一个3×3矩阵的例子。已知 C ＝ 。在计算行列式时需要用到子式（minor）和余子式（cofactor）。矩阵元素 的子式是删除第 i 排和第 j 列后得到的矩阵行列式。例如，在矩阵 C 中，顺着第一排展开， a ₁₁ ＝2的子式是行列式 ＝7·6－8·9＝－30， a ₁₂ ＝1的子式是行列式＝3·6－4·9＝－18， a ₁₃ ＝5的子式是行列式 ＝3·8－4·7＝－4。 a _ij 的余子式是 ×子式，因而余子式是带符号的子式，有时与子式相同，有时与子式相反。本例 a ₁₁ 、 a ₁₂ 和 a ₁₃ 的余子式为－30、18和－4。矩阵行列式为各个元素与对应的余子式乘积之和。针对本例， 2·（－30）＋1·（18）＋5·（－4）＝－62。如果按照第一列展开， ＝2的子式是－30， ＝3的子式是行列式 ＝1·6－8·5＝－34， a ₃₁ ＝4的子式是行列式 ＝1·9－7·5＝－26。它们对应的余子式为－30、34和－26， 。不论按排还是按列拓展余子式，矩阵的行列式不变。由此可以看出，行列式是方阵的重要性质。R计算行列式的函数是det，如det（C）＝－62。

行列式在多元（multivariate）统计分析中的重要性表现在两个相关的方面（Pituch＆ Stevens，2016，p.52）。其一，协方差矩阵（covariance matrix）的行列式是多元变量的广义方差（generalized variance），表示在排除一组变量共同方差之后还剩下多少变异性（variability）。其二，行列式是一组变量的方差测量，几个多元统计量（譬如多元方差分析中的Wilks’ Λ ）的计算都需要行列式。

逆矩阵 （ inverse of a matrix ）。方阵 X 的逆矩阵 X ^－1 满足以下条件： XX ^－1 ＝ X ^－1 X ＝ I （单位阵）。这里以 C ＝ 为例介绍逆矩阵算法。

先计算矩阵的子式和余子式矩阵。通过计算得到矩阵 C 的子式矩阵 ，对应的余子式为 。然后，计算余子式的转置矩阵，转置后的矩阵称作伴随矩阵（adjoint matrix）。本例伴随矩阵为 。最后，将伴随矩阵中的每个元素除以行列式，得到矩阵的逆矩阵。本例矩阵 C 的行列式为－62，因此矩阵 C 的逆矩阵为：

检验矩阵 C ^－1 是矩阵 C 的逆矩阵：

R计算逆矩阵的函数是solve，如solve（C）。矩阵没有除法，但是一个矩阵乘以另一个矩阵的逆矩阵类似于矩阵相除。如果一个方阵 X 的逆矩阵存在，那么 X 就是非奇异矩阵（nonsingular matrix），否则 X 就是奇异矩阵（singular matrix）。

秩 （ rank ）。 X 是 m × n 矩阵，有长度为 m 的 n 列（ x ₁ ， x ₂ ，…， x _n ）向量，长度为 n 的 m 排向量。如果一组向量（ x ₁ ， x ₂ ，…， x _k ）线性组合（linear combinations）不为0，则该组向量是线性独立的（linearly independent），否则有线性相依性。矩阵的秩是行列式不为0的 X 最大子方阵（submatrix）的阶（order）。 X 的列秩（column rank）是该矩阵线性独立列的最大数，记作C（A）；排秩（row rank）是该矩阵线性独立排的最大数，记作R（ A ）。矩阵的秩是最大子方阵数，因而列秩和排秩是相同的，通常只说矩阵秩（ ρ （ X ）），无需具体地说排秩或列秩。矩阵 A 的秩总是小于或等于 m 和 n 中的小值，即 ρ （ X ）≤ min （ m ， n ）。当 m × n 矩阵 X 的秩 ρ （ X ）＝ min （ m ， n ）时，矩阵为满秩（full rank）。判断一个 n × n 方阵是否为满秩，最简单的方法是看行列式。只有当方阵为满秩时，方阵的行列式才不会为0。

这里举几个求矩阵秩的例子。

令 X ＝ 。 X 是2×3矩阵，因此 ρ （ X ）≤2， ρ （ X ）要么为1，要么为2。如果 ρ （ X ）＝2，则有常数 a ₁ 和 a ₂ 使得 a ₁ （1，3，5）＋ a ₂ （2，4，7）＝0，那么 a ₁ ＋2 a ₂ ＝0，3 a ₁ ＋4 a ₂ ＝0，5 a ₁ ＋7 a ₂ ＝0。该组方程的解为 a ₁ ＝ a ₂ ＝0，矩阵 X 的两排向量是独立的，因此秩为2，即 ρ （ X ）＝2。

令 X ＝ 。 X 是2×2方阵，因此 ρ （ X ）≤2。这里，行列式为0，说明列向量线性相依，因此秩为1，即 ρ （ X ）＝1。

令 X ＝ 。 X 是3×4矩阵，因此 ρ （ X ）≤3。很容易看出， X 各排是线性独立的，否则一排向量元素会是另一排向量元素的乘数（乘数可以是整数或分数）。因此，矩阵 X 秩为3，即 ρ （ X ）＝3。

在R中安装数据包fBasics，调用函数rk（x，method＝'qr' ），可以计算任何矩阵的秩。也可以使用R内置函数qr（x）计算矩阵秩。矩阵秩在多元统计中发挥着重要作用。譬如，在多元线性回归模型中，如果矩阵 X 为非满秩，则表明变量之间有完全共线性存在。

特征值和特征向量 （ eigenvalue and eigenvector ）。特征值、特征向量与行列式一样都是方阵的基本属性。令 A 是 n × n 方阵，则 A 的特征值为矩阵 A 特征方程的特征根（roots of the characteristic equation）。特征值 λ 是标量。矩阵 A 特征方程为｜ A － λ I _n ｜＝0。即， A － λ I _n 是奇异矩阵，其行列式值为0。满足 Ax ＝ λ x 的向量 x 称作特征向量。

以矩阵E＜－ matrix（c（0，2，3，1），ncol＝2）为例求特征根。 － ＝0。根据行列式性质， 。解方程求得： ＝3， ＝－2。当 ＝3时， ＝3 ，则 ＝ 。当 ＝－2时， ＝－2 ，则 。R 计算特征根和特征向量的函数为eigen。对于本例，执行R命令eigen（E）后的结果为：

R输出的特征值是正规解（normalized solutions）。正规解满足 x′x ＝1。譬如， ＝ _时 ， ＝1，经计算得到 。正规化向量乘以－1不会改变向量的属性，因此 ≈－0.707 106 8。

特征值和特征向量有以下特点：①矩阵特征值之和等于矩阵迹；②矩阵特征值乘积等于矩阵行列式；③对于一个非奇异矩阵，特征值数等于矩阵的秩（Hagle，1995，pp.90－91）。譬如，上例矩阵 E 特征值之和等于1，矩阵迹也为1；矩阵 E 特征值乘积等于－6，行列式也为－6；矩阵 E 特征值数为2，矩阵的秩也为2。

特征值和特征向量用于主成分分析（principal component analysis）、因子分析（factor analysis）、回归诊断（regression diagnosis）和典型相关分析（canonical correlation analysis）等多元统计中。譬如，在多元回归模型中，特征值的大小表示变量的关联度。当多元回归模型中矩阵 X′X 是非奇异矩阵，但是至少有一个很小的特征值时，变量之间接近多元共线性（multicollinearity）（Gruber，2014，p.143）。

思考与练习

1.简要说明R向量、数据框和矩阵之间的关系。

2.列举一些常用的特殊矩阵，简要说明它们各自的特点。

3.R变量x储存数据集（6，31，38，24，30，31，34，29，40，29，32，35，28，25，39，18，34，30，35，32），利用R自带函数sum和mean计算x数值之和与平均数。

4.针对第3题，如果不使用R自带函数mean计算x数值的平均数，计算平均数的R命令是什么？

5.针对第3题，编写R命令得到x＞30的所有数值，并计算这些数值的总个数及数值之和。

6.有以下两个矩阵：

利用R命令计算矩阵 X 和 Y 的点乘积和克罗内克积。

7.有两组数据，第一组（G1）在因变量上的测量数值为：33，41，53，31，39，41，46，38，56，39，43，48，37，32，54，22，47，40，48，43；第二组（G2）在同一个因变量上的测量数值为：66，33，61，65，45，20，50，39，53，48，52，48，56，42，37，39，28，36，39，43。利用这两组数据创建一个数据框，组别变量名为G，组别变量水平为G1和G2，因变量名称为DV，数据框名称为Data。

8.编写R命令计算第7题数据框Data中每个组的平均数。

9.已知方阵 A ＝ ，编写R命令计算方阵 A 的迹、行列式、逆矩阵和特征值。

10.已知以下三个数值向量：

Class1＜－c（27，17，12，23，13，19，25，16，11，17，34）；

Class2＜－c（19，18，25，36，17，28，20，14，15，10，22）；

Class3＜－c（17，14，25，16，31，18，22，33，14，18，20）。

编写R命令创建一个数据列表，名称为Class，利用R函数lapply计算每个列表数据的平均数。

11.已知一元二次方程 ＋ bx ＋ c ＝0 两个实根的计算公式为： x ＝ ，其中判别式 －4 ac ≥0 。要求利用这个公式编写R 函数quadratic ，并利用该函数计算方程 －6 x －20＝0的两个实根。 YzdGu4nWkm75hMxjYu1UsDMq2gPh2Lev7N676dFrOOUYnXb1RFlUMyuUO3jsDkwH