如何在NumPy/SciPy中高效对角化三角矩阵？

核心思路 #

三角矩阵的特征值可直接提取对角元素，这一步完全无需调用通用特征值求解函数；特征向量则可利用三角矩阵的结构直接递推计算，比numpy.linalg.solve或通用eig函数更高效、精准。

1. 特征值的快速获取 #

直接提取矩阵对角元素即可，这是O(n)操作，速度远快于通用特征值算法：

import numpy as np

# 假设tri_mat是三角矩阵（上/下三角均可）
eigenvalues = np.diag(tri_mat)

2. 特征向量的高效递推计算 #

利用三角矩阵的稀疏结构，无需调用通用线性方程组求解器，直接逐行递推得到特征向量：

上三角矩阵版本

对于上三角矩阵，特征值λ对应第k个对角元素时，特征向量的第k个分量设为1，从第k-1行往前递推求解其余分量：

def upper_triangular_eigvecs(mat):
    n = mat.shape[0]
    eigenvalues = np.diag(mat)
    eigenvectors = np.zeros((n, n), dtype=mat.dtype)
    
    for k in range(n):
        lam = eigenvalues[k]
        vec = np.zeros(n, dtype=mat.dtype)
        vec[k] = 1.0
        # 从k-1行反向递推到第0行
        for i in range(k-1, -1, -1):
            denom = mat[i, i] - lam
            # 单重特征值下分母非零，重特征值需额外处理
            if not np.isclose(denom, 0):
                vec[i] = -np.dot(mat[i, i+1:n], vec[i+1:n]) / denom
        eigenvectors[:, k] = vec
    return eigenvalues, eigenvectors

下三角矩阵版本

类似地，下三角矩阵从第k+1行正向递推：

def lower_triangular_eigvecs(mat):
    n = mat.shape[0]
    eigenvalues = np.diag(mat)
    eigenvectors = np.zeros((n, n), dtype=mat.dtype)
    
    for k in range(n):
        lam = eigenvalues[k]
        vec = np.zeros(n, dtype=mat.dtype)
        vec[k] = 1.0
        # 从k+1行正向递推到第n-1行
        for i in range(k+1, n):
            denom = mat[i, i] - lam
            if not np.isclose(denom, 0):
                vec[i] = -np.dot(mat[i, 0:i], vec[0:i]) / denom
        eigenvectors[:, k] = vec
    return eigenvalues, eigenvectors

3. 批量矩阵处理优化 #

针对200个200阶三角矩阵的场景，直接循环调用上述函数即可，单进程就能高效完成：

# 假设batch_mats是形状为(200, 200, 200)的3D数组，每个元素是三角矩阵
batch_eigenvalues = []
batch_eigenvectors = []

for mat in batch_mats:
    # 根据矩阵是上/下三角选择对应函数
    if np.allclose(mat, np.triu(mat)):
        vals, vecs = upper_triangular_eigvecs(mat)
    else:
        vals, vecs = lower_triangular_eigvecs(mat)
    batch_eigenvalues.append(vals)
    batch_eigenvectors.append(vecs)

# 转换为numpy数组方便后续处理
batch_eigenvalues = np.array(batch_eigenvalues)
batch_eigenvectors = np.array(batch_eigenvectors)

4. 与numpy.linalg.eig的对比 #

numpy的eig函数并未公开说明对三角矩阵的优化逻辑，其内部仍会执行通用的QR迭代流程，远不如上述利用三角结构的方法高效。同时，手动递推避免了通用算法中的额外精度损耗，结果更精准。

重特征值说明 #

若三角矩阵存在重特征值，上述方法仅适用于矩阵可对角化的场景（即重特征值对应的Jordan块为1x1）。若需处理不可对角化的三角矩阵，需额外构造Jordan基，但此类场景下矩阵本身无法对角化，需调整需求。

内容的提问来源于stack exchange，提问作者Nim