自定义行约简：将矩阵前N/2列转为N/2阶单位矩阵

问题描述 #

现有一个k行p列的Vandermonde矩阵G，所有元素均为模素数p的结果（p略大于k）。需要通过行初等变换，将矩阵前m行的前m列转化为m×m单位矩阵（m = k//2），要求：

• 不得使用sympy.rref工具（该工具会自动进行行交换及全矩阵化简，不符合需求）
• 仅完成前m列的化简即停止，其余行保持原样

生成该Vandermonde矩阵的代码如下：

import numpy as np

def ReedSolomon(k,p):
    rr = np.arange(k).reshape((-1, 1))
    cc = np.arange(p)
    return cc**rr % p

示例 #

调用ReedSolomon(4,5)得到原矩阵：

[[1, 1, 1, 1, 1],
 [0, 1, 2, 3, 4],
 [0, 1, 4, 4, 1],
 [0, 1, 3, 2, 4]]

期望化简后矩阵（前2行前2列为单位矩阵）：

[[1, 0, -1, -2, -3],
 [0, 1, 2, 3, 4],
 [0, 1, 4, 4, 1],
 [0, 1, 3, 2, 4]]

注：负数在模p下等价于正数（如-1 mod5=4），代码实现中可根据需求选择保留负数或转为正数。

实现思路 #

由于目标矩阵的前m行前m列是上三角矩阵（第一列仅第一行非零，对角线元素均与p互质），无需行交换即可完成化简：

1. 遍历每一个主元列（从0到m-1）
2. 将主元行的主元元素化为1（通过乘以模p逆元）
3. 消去主元列上方所有行的对应元素，最终得到单位矩阵子块

代码实现 #

import numpy as np

def ReedSolomon(k,p):
    rr = np.arange(k).reshape((-1, 1))
    cc = np.arange(p)
    return cc**rr % p

def partial_row_reduce(mat, p):
    # 复制原矩阵，避免修改输入
    mat_copy = mat.copy().astype(np.int64)
    m = mat_copy.shape[0] // 2  # 需要化简的前m行前m列
    for c in range(m):
        # 将主元元素转为1：计算模p逆元（费马小定理）
        pivot = mat_copy[c, c]
        if pivot != 1:
            inv_pivot = pow(pivot, p-2, p)
            mat_copy[c] = (mat_copy[c] * inv_pivot) % p
        # 消去当前列上方的所有非零元素
        for r in range(c):
            factor = mat_copy[r, c]
            if factor != 0:
                mat_copy[r] = (mat_copy[r] - factor * mat_copy[c]) % p
    return mat_copy

# 测试示例
if __name__ == "__main__":
    k = 4
    p = 5
    original_mat = ReedSolomon(k, p)
    print("原矩阵：")
    print(original_mat)
    reduced_mat = partial_row_reduce(original_mat, p)
    print("\n化简后矩阵：")
    print(reduced_mat)

代码说明 #

• partial_row_reduce函数接收原矩阵和素数p，返回化简后的矩阵
• 利用费马小定理计算模p逆元：对于素数p，任意与p互质的a，其逆元为a^(p-2) mod p
• 所有运算均在模p下进行，确保结果符合素数域要求
• 仅处理前m行的前m列，其余行保持原状态，避免过度化简

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