46、稀疏矩阵基本操作与特征值特征向量计算

稀疏矩阵基本操作与特征值特征向量计算

1. 稀疏矩阵的基本操作

1.1 ij 形式的稀疏矩阵

对于无权图的邻接矩阵，如果所有元素的值都相等，那么只需要指定行索引向量 i 和列索引向量 j 即可。若图是有向图，向量 i 、 j （必要时还有向量 s ，用于存储非零元素的值）的长度等于弧的数量 K ；若图是无向图，邻接矩阵 A 是对称的，理论上只需存储下三角部分（即 i > j 的非零元素），但为了方便某些操作（如遍历节点的邻居），通常会存储矩阵 A 所有非零元素的索引，此时向量 i 、 j 的长度是边数量的两倍，即 M = 2K 。

使用 ij 形式进行矩阵 - 向量乘法比使用完整的邻接矩阵更高效。设 x 是一个 N 维列向量，矩阵 - 向量乘积 y = Ax 可以按以下算法计算：

Algorithm 9 Matrix vector multiplication (ij-form)
Input: i, j, s, x
Output: y = Ax
1: y ←⃗0
2: for k = 0 to M −1 do
3:
    y[i[k]] ←y[i[k]] + s[k] ∗x[j[k]]
4: end for

该算法执行 M 次乘法和 N 次初始化操作（将向量 y 的所有分量初始化为零），因此所需操作数为 M + N ，时间复杂度为 $O(M)$。而使用完整邻接矩阵进行矩阵 - 向量乘法需要两层嵌套循环，时间复杂度为 $O(N^2)$。由于稀疏矩阵中 M ≪ N^2 ，所以 ij 形式的矩阵 - 向量乘法更高效。

1.2 压缩行存储（CRS）格式

另一种常见的稀疏矩阵存储方式是压缩行存储（CRS）格式。在 CRS 格式中，矩阵信息存储在三个向量 r 、 j 、 s 中。与 ij 形式不同，这三个向量的长度不同，其中 j 和 s 有 M 个分量， r 有 N + 1 个分量。向量 j 按行索引顺序存储矩阵每个非零元素的列索引，向量 s 存储矩阵的非零值，向量 r 中， r[i] 表示第 i 行第一个非零元素在 j 和 s 中的索引，且 r[N] = M 。这样，矩阵第 i 行的所有非零元素将存储在位置 k = r[i], r[i] + 1, r[i] + 2, ..., r[i + 1] - 1 中。

以下是使用 CRS 格式进行矩阵 - 向量乘法的算法：

Algorithm 10 Matrix vector multiplication (CRS-form)
Input: r, j, s, x
Output: y = Ax
1: y ←⃗0
2: for i = 0 to N −1 do
3:
    for k = r[i] to r[i + 1] −1 do
4:
        y[i] ←y[i] + s[k] ∗x[j[k]]
5:
    end for
6: end for

有时候需要在不同的稀疏矩阵表示形式之间进行转换。从 CRS 格式转换为 ij 格式很简单，只需按以下算法定义向量 i ：

Algorithm 11 From CRS to ij-format
Input: r
Output: i
1: for ℓ= 0 to N −1 do
2:
    for k = r[ℓ] to r[ℓ+ 1] −1 do
3:
        i[k] ←ℓ
4:
    end for
5: end for

从 ij 格式转换为 CRS 格式稍微复杂一些，因为图的边在 i 和 j 向量中没有特定的顺序。一种方法是先对边 (i, j) 按 i 的升序排序，然后扫描排序后的向量 i 来创建向量 r ，但有更简单高效的方法，如下所示：

Algorithm 12 From ij-to CRS format
Input: i, j
Output: ⃗r, jCRS
1: ⃗r ←⃗0, ⃗p ←⃗0 {Set ⃗r ∈NN+1 and ⃗p ∈NN to zero}
2: Count the number of elements for each row
3: for k = 0 to M −1 do
4:
    r[i[k] + 1] ←r[i[k] + 1] + 1
5: end for
6: Compute vector r
7: r[0] ←0
8: for l = 0 to N −1 do
9:
    r[l + 1] ←r[l] + r[l + 1]
10: end for
11: Compute vector jCRS
12: for k = 0 to M −1 do
13:
    ℓ←i[k]
14:
    jCRS[r[ℓ] + p[ℓ]] ←j[k]
15:
    p[ℓ] ←p[ℓ] + 1
16: end for

这个过程可以通过一个示例更好地理解。假设向量 i 的前 11 个元素中，标签 0 出现 4 次（ p0 = 4 ），标签 1 出现 3 次（ p1 = 3 ），标签 2 出现 2 次（ p2 = 2 ），标签 3 出现 2 次（ p3 = 2 ），且在算法第 7 - 10 行生成的向量 r = (0, 7, 14, 18, 21, ...) ，这意味着节点 0 有 r1 - r0 = 7 条边，节点 1 有 r2 - r1 = 7 条边，节点 2 有 r3 - r2 = 4 条边，节点 3 有 r4 - r3 = 3 条边。当读取第 12 条边时，如果 i[11] = 1 且 j[11] = 17 ，则将值 17 插入 jCRS[r[1] + p[1]] 位置，并将 p1 加 1。

2. 特征值和特征向量计算

2.1 矩阵的特征值和特征向量

设 A 是一个 N × N 的实方阵，我们希望找到一个非零向量 u ，使得矩阵 A 作用于 u 后方向不变，即满足 Au = λu ，其中 λ 是一个标量。如果 λ 和 u 存在，则分别称为矩阵 A 的特征值和右特征向量（简称特征向量）。该方程可以改写为 (λI - A)u = 0 ，其中 I 是 N × N 的单位矩阵。由于这是一个齐次线性方程组，当 det(λI - A) ≠ 0 时，只有零解 u = 0 ；若要得到非零解，则需满足 p(λ) ≡ det(λI - A) = 0 ，其中 p(λ) 称为矩阵 A 的特征多项式。矩阵 A 的特征值就是特征多项式 p(λ) 的根，根据代数基本定理， p(λ) 是一个 N 次多项式，有 N 个根（通常为复数）。找到特征值 λ1, ..., λN 后，可以通过求解线性齐次方程组 (λαI - A)uα = 0 （ α = 1, ..., N ）来计算特征向量。

矩阵 A 的所有特征值组成的集合 λ(A) = {λ1, ..., λN} 称为矩阵 A 的谱，所有特征值的最大绝对值称为矩阵的谱半径，记为 ρ(A) = max1≤α≤N |λα| 。

类似地，可以定义左特征向量 v ，满足 v†A = λv† ，其中 † 表示向量的伴随（对于列向量 x ，其伴随 x† 是一个行向量，且 x† 的第 i 个元素是 x 的第 i 个元素的复共轭）。可以证明，对应不同特征值的左、右特征向量是正交的，即若 λα ≠ λβ ，则 v†αuβ = 0 。

2.2 可对角化矩阵

实对称矩阵以及更一般的 Hermitian 矩阵（满足 A = A† 的方阵，其中 A† 是 A 的转置的复共轭）在特征值和特征向量方面有特殊性质。实对称矩阵或 Hermitian 矩阵的特征值是实数。以 Au = λu 为例，取其伴随得到 u†A† = λu† ，将第一个式子左乘 u† ，第二个式子右乘 u 后相减，可得 u†(A - A†)u = (λ - λ)u†u 。由于 u†u = ∥u∥2 > 0 ，当 A = A† 时，有 λ = λ ，即 λ 是实数。同时，实对称矩阵或 Hermitian 矩阵的左、右特征向量相同，且特征向量也是实数，因为它们是实线性齐次方程组的解。

通常会将特征向量归一化，使其欧几里得范数等于 1，即 u⊤αuα = 1 （ α = 1, ..., N ）。如果所有特征值都不同（即 λα ≠ λβ ， α ≠ β ， α, β = 1, ..., N ），则有 u⊤αuβ = δαβ ，用矩阵形式表示为 U⊤U = I ，其中 U = (u1, ..., uN) 是由归一化特征向量作为列向量组成的方阵。这表明特征向量是线性无关的，且 U 的转置等于其逆，即 U−1 = U⊤ 。

对于非对称矩阵，其特征值不一定是实数，但如果所有特征值都不同，则对应的 N 个特征向量是线性无关的，它们构成 CN 的一个基。以下定理说明了对应不同特征值的特征向量是线性无关的：
定理 A.1 ：设 λ1, λ2, ..., λk 是矩阵 A ∈ CN×N 的 k ≤ N 个不同特征值， u1, ..., uk 是对应的特征向量，则这 k 个特征向量线性无关。

如果矩阵 A 的特征向量构成 CN 的一个基，则 det(U) ≠ 0 ，可以将矩阵 A 对角化，即 Λ = U−1AU ，其中 Λ = diag(λ1, ..., λN) 是包含特征值的对角矩阵。如果特征值不全不同，矩阵是否可对角化取决于矩阵的结构。设不同的特征值为 λ1, ..., λm ，特征值 λi 的重数为 mi ，且 ∑mi = N ，矩阵可对角化的充要条件是对于每个特征值 λi 都有恰好 mi 个线性无关的特征向量。如果不满足该条件，矩阵不能通过相似变换化为对角矩阵，但可以化为 Jordan 标准形。

2.3 幂法

幂法和 Lanczos 法是计算矩阵最大特征值及其对应特征向量的两种常用方法。虽然基于 Lanczos 法的算法可能更适合计算稀疏矩阵的特征值和特征向量，但这里主要介绍幂法。

设 A 是一个 N × N 的实方阵，其特征值按绝对值降序排列为 λ1, ..., λN ，且假设 |λ1| > |λ2| ，即 |λ1| > |λ2| ≥ |λ3| ≥ ... ≥ |λN| 。这意味着 λ1 是单重且实数的（可以用反证法证明，如果 λ1 有非零虚部，由于矩阵的特征多项式系数为实数， λ1 的共轭也是特征值，这会导致两个最大特征值的绝对值相等，与假设矛盾），对应的特征向量 u1 也是实数。

考虑向量序列 xn+1 = Axn （ n = 0, 1, ... ），其中 x0 是 RN 中的任意非零向量，例如 x0 = (1, 1, ..., 1)⊤ 。可以证明：
$\lim_{n \to \infty} \frac{x_n}{|x_n|} = u_1$
即归一化后的向量 xn/∥xn∥ 收敛到与主特征值 λ1 对应的归一化特征向量 u1 ，且当 n → ∞ 时， xn+1 与 xn 成比例，比例常数为 λ1 。

由于矩阵 A 的特征向量 (u1, u2, ..., uN) 构成 RN 的一个基，任意向量 x0 可以表示为这些特征向量的线性组合：
$x_0 = \sum_{i = 1}^{N} \alpha_i u_i$
$x_1 = Ax_0 = A\sum_{i = 1}^{N} \alpha_i u_i = \sum_{i = 1}^{N} \alpha_i \lambda_i u_i$
$x_2 = Ax_1 = A\sum_{i = 1}^{N} \alpha_i \lambda_i u_i = \sum_{i = 1}^{N} \alpha_i \lambda_i^2 u_i$
一般地，有：
$x_n = \lambda_1^n (\alpha_1 u_1 + \sum_{i = 2}^{N} \alpha_i (\frac{\lambda_i}{\lambda_1})^n u_i)$
如果 α1 ≠ 0 ，由于 |λ1| > |λi| （ i ≠ 1 ），则当 n → ∞ 时，$(\frac{\lambda_i}{\lambda_1})^n → 0$，这表明 xn 会与第一个特征向量 u1 成比例。

考虑 xn+1 和 xn 的第 j 个分量的比值：
$\frac{x_{n + 1, j}}{x_{n, j}} = \lambda_1 \frac{\alpha_1 u_{1, j} + \sum_{i = 2}^{N} \alpha_i (\frac{\lambda_i}{\lambda_1})^{n + 1} u_{i, j}}{\alpha_1 u_{1, j} + \sum_{i = 2}^{N} \alpha_i (\frac{\lambda_i}{\lambda_1})^n u_{i, j}}$
当 α1 u1,j ≠ 0 时，随着 n → ∞ ，该比值会收敛到 λ1 。实际上，该比值渐近地以公比为 $\frac{\lambda_2}{\lambda_1}$ 的几何序列收敛到 λ1 ，即：
$\frac{x_{n + 1, j}}{x_{n, j}} = \lambda_1 + O((\frac{\lambda_2}{\lambda_1})^n)$
这意味着每次迭代时， λ1 的估计误差会以常数因子 $|\frac{\lambda_2}{\lambda_1}|$ 减小，因此 $\frac{\lambda_2}{\lambda_1}$ 越小，比值 $\frac{x_{n + 1, j}}{x_{n, j}}$ 收敛到 λ1 的速度越快。

为了更好地利用 xn 和 xn+1 的信息，引入 Rayleigh 商：
$\sigma_n = \frac{\sum_{j} x_{n, j} x_{n + 1, j}}{\sum_{j} x_{n, j} x_{n, j}} = \frac{x_n^⊤ Ax_n}{x_n^⊤ x_n}$
Rayleigh 商有一个重要性质：对于任意近似特征向量 x ，使残差 η = Ax - ˜λx 最小的 ˜λ 就是 Rayleigh 商。可以证明，一般情况下 $\sigma_n = \lambda_1 + O((\frac{\lambda_2}{\lambda_1})^n)$，当矩阵 A 是对称矩阵时，有 $\sigma_n = \lambda_1 + O((\frac{\lambda_2}{\lambda_1})^{2n})$，即对称矩阵的 Rayleigh 商收敛速度更快。

然而，简单地迭代 xn+1 = Axn 在大多数情况下会出现问题。如果 |λ1| > 1 ， xn 会指数级发散，可能导致机器溢出；如果 |λ1| < 1 ， xn 会指数级收敛到 0，可能导致机器下溢。为避免这个问题，在每次迭代时对 xn 进行重新归一化，定义新的向量序列：
$y_n = \frac{x_n}{|x_n|}$
$x_{n + 1} = Ay_n$

最后，需要确定何时停止迭代，认为 σn 对 λ1 的近似是满意的。一种实用的停止准则是要求 σn 对 λ1 的近似具有较小的相对残差，即 ∥σn - σn−1∥ ≤ ε∥σn∥ ，其中 ε 是预先指定的容差。但这个准则在收敛速度很慢时可能有风险，更好的停止准则基于残差：如果矩阵 A 有一组线性无关的特征向量，且 ∥Ayn - σnyn∥2 ≤ ε∥Ayn∥2 ，则有 $\min_{j} |1 - \frac{\sigma_n}{\lambda_j}| ≤ ε\mu_2(U)$，其中 $\mu_2(U) = |U|_2 |U^{-1}|_2$ 是特征向量矩阵 U 在欧几里得范数下的条件数。当矩阵 A 是对称矩阵时， U 是酉矩阵，有 $|U|_2 = |U^{-1}|_2 = 1$，则 $|1 - \frac{\sigma_n}{\lambda_1}| ≤ ε$。

幂法的伪代码如下：

Algorithm 13 Power method
Input: A, ϵ
Output: u1, λ1
1: x0 ←(1, ..., 1)⊤∈RN
2: n ←0
3: yn ←xn/∥xn∥
4: xn+1 ←Ayn
5: σn ←y⊤n xn+1
6: if ∥xn+1 −σnyn∥> ε∥xn+1∥then
7:
    n ←n + 1
8:
    goto 3
9: else
10:
    λ1 ←σn
11:
    u1 ←yn
12: end if

幂法的实现需要存储两个向量 xn 和 yn 以及稀疏矩阵 A ，每次迭代执行一次矩阵 - 向量乘法（时间复杂度为 $O(K)$）和两次范数计算（时间复杂度为 $O(N)$），因此算法的时间复杂度为 $O(K × n_ε)$，其中 n_ε 是在给定误差 ε 下近似 λ1 和 u1 所需的（未知）步数。

综上所述，稀疏矩阵的存储和运算以及特征值和特征向量的计算在数值分析和科学工程中有重要应用。不同的稀疏矩阵存储格式和算法在效率和适用性上有所不同，需要根据具体问题选择合适的方法。幂法是一种简单有效的计算矩阵主特征值和对应特征向量的方法，但在实际应用中需要注意数值稳定性和收敛速度等问题。

3. 方法总结与对比

3.1 稀疏矩阵存储格式对比

存储格式	存储向量	向量长度	适用场景	矩阵 - 向量乘法复杂度	转换复杂度
ij 形式	i , j , s	i , j , s 长度为 M	图的邻接矩阵表示，操作简单，适合快速实现	$O(M)$	从 CRS 转 ij 简单，从 ij 转 CRS 较复杂
压缩行存储（CRS）格式	r , j , s	j , s 长度为 M ， r 长度为 N + 1	适合大规模稀疏矩阵，便于按行访问元素	$O(M)$	从 CRS 转 ij 简单，从 ij 转 CRS 有高效算法

3.2 特征值计算方法对比

方法	适用场景	收敛速度	复杂度	稳定性
幂法	计算矩阵的主特征值和对应特征向量	取决于 $\vert\frac{\lambda_2}{\lambda_1}\vert$，该值越小收敛越快；对称矩阵收敛更快	$O(K × n_ε)$， n_ε 为迭代步数	需注意数值溢出和下溢问题，可通过归一化解决
Lanczos 法	适合计算稀疏矩阵的特征值和特征向量	-	-	未详细介绍，可能更高效但超出本文范围

3.3 流程图总结

graph TD;
    classDef startend fill:#F5EBFF,stroke:#BE8FED,stroke-width:2px;
    classDef process fill:#E5F6FF,stroke:#73A6FF,stroke-width:2px;
    classDef decision fill:#FFF6CC,stroke:#FFBC52,stroke-width:2px;

    A([开始]):::startend --> B(选择稀疏矩阵存储格式):::process;
    B --> C{是否为 ij 形式}:::decision;
    C -- 是 --> D(进行 ij 形式矩阵 - 向量乘法):::process;
    C -- 否 --> E(进行 CRS 形式矩阵 - 向量乘法):::process;
    D --> F{是否需要转换格式}:::decision;
    E --> F;
    F -- 是 --> G{从 ij 转 CRS 还是从 CRS 转 ij}:::decision;
    F -- 否 --> H(是否计算特征值和特征向量):::decision;
    G -- 从 ij 转 CRS --> I(使用 Algorithm 12 转换):::process;
    G -- 从 CRS 转 ij --> J(使用 Algorithm 11 转换):::process;
    I --> H;
    J --> H;
    H -- 是 --> K(选择特征值计算方法):::process;
    H -- 否 --> L([结束]):::startend;
    K --> M{是否使用幂法}:::decision;
    M -- 是 --> N(使用幂法迭代计算):::process;
    M -- 否 --> O(使用其他方法，如 Lanczos 法):::process;
    N --> P{是否满足停止准则}:::decision;
    O --> P;
    P -- 是 --> Q(输出特征值和特征向量):::process;
    P -- 否 --> N;
    Q --> L;

4. 实际应用中的考虑因素

4.1 稀疏矩阵存储格式选择

在实际应用中，选择稀疏矩阵的存储格式需要考虑以下因素：
- 数据特点 ：如果图的边没有特定顺序，ij 形式可能更方便初始存储；如果需要频繁按行访问元素，CRS 格式更合适。
- 操作需求 ：如果只是简单的矩阵 - 向量乘法，两种格式都可以；如果需要进行格式转换，需要考虑转换的复杂度。
- 内存使用 ：CRS 格式可能需要额外的 r 向量，但在某些情况下可以更高效地利用内存。

4.2 特征值计算方法选择

计算矩阵的特征值和特征向量时，选择合适的方法至关重要：
- 矩阵规模 ：对于大规模稀疏矩阵，Lanczos 法可能更高效，但实现复杂；幂法简单易实现，适合小规模矩阵或只需要主特征值的情况。
- 收敛速度 ：幂法的收敛速度取决于 $\vert\frac{\lambda_2}{\lambda_1}\vert$，如果该值较大，收敛可能较慢，需要考虑其他方法。
- 稳定性 ：幂法在迭代过程中可能出现数值溢出或下溢问题，需要进行归一化处理；其他方法可能有更好的稳定性，但需要进一步研究。

5. 代码实现示例（Python 伪代码）

5.1 ij 形式矩阵 - 向量乘法

def matrix_vector_multiplication_ij(i, j, s, x):
    N = len(x)
    M = len(i)
    y = [0] * N
    for k in range(M):
        y[i[k]] += s[k] * x[j[k]]
    return y

5.2 CRS 形式矩阵 - 向量乘法

def matrix_vector_multiplication_crs(r, j, s, x):
    N = len(x)
    y = [0] * N
    for i in range(N):
        for k in range(r[i], r[i + 1]):
            y[i] += s[k] * x[j[k]]
    return y

5.3 幂法实现

import numpy as np

def power_method(A, epsilon):
    N = A.shape[0]
    x0 = np.ones(N)
    n = 0
    while True:
        xn = x0 / np.linalg.norm(x0)
        xn_plus_1 = A @ xn
        sigma_n = np.dot(xn, xn_plus_1)
        if np.linalg.norm(xn_plus_1 - sigma_n * xn) > epsilon * np.linalg.norm(xn_plus_1):
            x0 = xn_plus_1
            n += 1
        else:
            lambda_1 = sigma_n
            u1 = xn
            break
    return u1, lambda_1

6. 总结与展望

本文详细介绍了稀疏矩阵的基本操作（包括 ij 形式和压缩行存储格式）以及特征值和特征向量的计算方法（主要是幂法）。稀疏矩阵的存储格式和特征值计算方法在数值分析、科学工程等领域有广泛应用，不同的方法在效率、复杂度和稳定性上各有优劣。

在未来的研究中，可以进一步探索以下方向：
- 更高效的稀疏矩阵存储格式和算法 ：随着数据规模的不断增大，需要开发更高效的存储和运算方法，以提高计算效率和减少内存使用。
- 特征值计算方法的改进 ：研究如何提高幂法等算法的收敛速度和稳定性，或者开发新的特征值计算方法，以满足不同场景的需求。
- 实际应用案例研究 ：将这些方法应用到具体的科学和工程问题中，验证其有效性和实用性，并根据实际需求进行优化。

通过不断的研究和实践，我们可以更好地利用稀疏矩阵和特征值计算方法，解决各种复杂的实际问题。