基于Cholesky分解的对称带状线性方程组前向替换问题排查

问题根源：索引计算逻辑错误 #

根据你定义的对称带状存储规则A[i,j] = Ab[i-j,j]（j ≤ i < min(n,j+1+p)），下三角Cholesky因子L[i,j]（i≥j）在带状存储矩阵中的行索引应为i-j，而非代码中错误的j+i。这个索引错误导致当j=1、i=2时（对应测试案例中L[2,1]的位置），程序取到了完全无关的元素，最终x[2]未被正确更新，保留了初始的b[2]=6。

修正后的前向替换代码 #

for j in range(n):
    c = min(j + 1 + p, n)
    # 先处理当前行的对角元（带状存储的第0行对应对角元素）
    x[j] = x[j] / L[0, j]
    # 遍历当前列下方的带状非零元素
    for i in range(j + 1, c):
        # 正确索引：i-j是带状存储的行号，j是列号
        x[i] = x[i] - L[i - j, j] * x[j]

测试案例验证细节 #

针对你的测试案例（p=1）：

• 非带状存储的Cholesky因子L为：

  [[3.46410162, 0, 0],
   [1.15470054, 0.81649658, 0],
   [0, -1.22474487, 2.44948974]]
  

• 对应带状存储的L（2行3列）为：

  [[3.46410162, 0.81649658, 2.44948974],  # 对角元行（i-j=0）
   [1.15470054, -1.22474487, 0]]          # 次对角元行（i-j=1）
  

修正索引后，程序会正确调用L[1,1]（即-1.2247）更新x[2]，最终得到预期的全1向量。

反向替换实现补充 #

反向替换用于求解Lᵀy = x，由于Lᵀ是上三角带状矩阵（上带宽等于原矩阵下带宽p），实现时需从后往前遍历：

for j in range(n-1, -1, -1):
    c = max(j - p, 0)
    # 先减去上方带状范围内的已计算值
    for i in range(c, j):
        # Lᵀ[i,j] = L[j,i]，对应带状存储的索引为j-i行、i列
        x[j] = x[j] - L[j - i, i] * x[i]
    # 处理对角元
    x[j] = x[j] / L[0, j]

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