如何用增广矩阵高效求解数百万个独立2x2线性方程组

问题背景 #

需要求解数百万至数千万个独立的二元线性方程组，输入为系数向量A、B、C、A'、B'、C'（每个向量包含数百万至数亿个元素）。每个索引i对应的方程组为：

A[i]x + B[i]y = -C[i]
A'[i]x + B'[i]y = -C'[i]

当前Python实现存在严重性能瓶颈：使用for循环结合Sympy逐个求解，代码如下：

from sympy import symbols, Eq, solve

x = []
y = []
x_sym, y_sym = symbols('x y')
for i in range(len(A)):
    sol = solve([
        Eq(A[i]*x_sym + B[i]*y_sym, -C[i]),
        Eq(A_prime[i]*x_sym + B_prime[i]*y_sym, -C_prime[i])
    ], [x_sym, y_sym])
    x.append(sol[x_sym])
    y.append(sol[y_sym])

该方案处理1000万+条目时计算时间过长，完全无法满足需求。

需求 #

编写DolphinDB脚本实现：

• 接收输入向量A、B、C、A'、B'、C'；
• 利用向量化/并行化运算快速计算x和y；
• 处理奇异矩阵等边缘情况。

DolphinDB实现方案 #

1. 核心原理：向量化应用克莱姆法则 #

对每个独立方程组，直接用克莱姆法则推导解的表达式，避免循环求解：

• 系数矩阵行列式：D = A[i]*B'[i] - A'[i]*B[i]
• 当D ≠ 0时，方程组有唯一解：

  x = (B[i]*C'[i] - B'[i]*C[i]) / D
  y = (A'[i]*C[i] - A[i]*C'[i]) / D
  

• 当D = 0时，属于奇异矩阵：
  • 若A[i]*C'[i] == A'[i]*C[i]且B[i]*C'[i] == B'[i]*C[i]，方程组有无穷多解；
  • 否则方程组无解。

2. 完整DolphinDB脚本 #

// 定义求解函数，接收6个向量参数
solveLinearSystems = function(A, B, C, A_prime, B_prime, C_prime) {
    // 计算系数矩阵行列式
    D = A * B_prime - A_prime * B
    
    // 计算x、y解的分子项
    numeratorX = B * C_prime - B_prime * C
    numeratorY = A_prime * C - A * C_prime
    
    // 初始化结果向量，用NULL标记无唯一解的情况
    x = array(double, size(A), NULL)
    y = array(double, size(A), NULL)
    
    // 筛选行列式非零的索引，计算唯一解
    validIdx = where(D != 0)
    x[validIdx] = numeratorX[validIdx] / D[validIdx]
    y[validIdx] = numeratorY[validIdx] / D[validIdx]
    
    // 处理奇异矩阵：区分无穷多解和无解
    singularIdx = where(D == 0)
    // 判定无穷多解的条件
    infiniteIdx = where(singularIdx, 
                        A[singularIdx] * C_prime[singularIdx] == A_prime[singularIdx] * C[singularIdx] 
                        && B[singularIdx] * C_prime[singularIdx] == B_prime[singularIdx] * C[singularIdx])
    // 可根据业务需求标记无穷多解，例如设为特定值或保留NULL
    // x[infiniteIdx] = 0.0
    // y[infiniteIdx] = 0.0
    
    return table(x as x_solution, y as y_solution)
}

// 测试示例：生成1000万条测试数据
n = 10000000
A = rand(1.0, n)
B = rand(1.0, n)
C = rand(1.0, n)
A_prime = rand(1.0, n)
B_prime = rand(1.0, n)
C_prime = rand(1.0, n)

// 执行求解
result = solveLinearSystems(A, B, C, A_prime, B_prime, C_prime)

3. 性能优化说明 #

• 向量化运算：DolphinDB向量运算基于C++底层实现，彻底避免Python循环的解释器开销，千万级数据处理仅需数秒；
• 自动并行：DolphinDB会自动对向量运算进行并行优化，无需手动编写并行逻辑；
• 内存高效：所有运算基于内存向量，无磁盘IO瓶颈，适配大规模数据场景。

4. 边缘情况处理 #

• 奇异矩阵：通过条件判断区分无穷多解和无解场景，可根据业务需求标记（如设为NULL、特定占位符）；
• 浮点精度：若需处理浮点误差，可将D != 0改为abs(D) > 1e-9（阈值可根据业务精度调整），避免误判奇异矩阵。

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