首先，咱们抓住问题的核心：你需要计算的x(1)本质上是单位向量e₁与A⁻¹b的点积，也就是x(1) = e₁ᵀA⁻¹b（其中e₁是第一个元素为1、其余为0的n维向量）。基于这个性质，我们可以用「预计算+快速点积」的思路，既保证精度，又大幅提升处理大量b时的效率。

具体步骤 #

• 预计算一次向量y：
  求解方程组Aᵀy = e₁（如果你的A是对称矩阵——毕竟你提到\调用LDL求解器，这类求解器通常针对对称矩阵，那Aᵀ=A，直接求解Ay = e₁即可）。在代码里可以这么写：

  n = size(A, 1);
  e1 = zeros(n, 1);
  e1(1) = 1;
  y = A \ e1;  % 或者A' \ e1，如果A非对称
  

  这一步只需要执行一次，因为A是固定不变的。而且这里用的是和A\b完全一致的稳定LDL求解器，不会有显式逆矩阵带来的精度损失。

• 对每个b快速计算x(1)：
  有了预计算好的y之后，每个b对应的x(1)就是y和b的点积：

  x1 = dot(y, b);  % 或者y' * b，注意维度匹配
  

  这一步的时间复杂度是O(n)，比每次调用A\b（即使利用缓存分解，也需要O(n²)的回代时间）高效得多，当你需要处理成百上千个b时，这个效率提升会非常显著。

为什么这个方案可行？ #

1. 精度保障：
   求解Ay=e1用的是和你平时求解Ax=b一样的LDL求解器，没有显式计算A的逆矩阵——而你已经提到显式逆矩阵会因为A的极差条件数损失精度，这个方法完美避开了这个问题，得到的x(1)精度和直接求解Ax=b后取第一个分量完全一致。

2. 效率提升：
   只需要对A做一次分解（第一次调用A\e1时完成，很多数值计算环境比如MATLAB会自动缓存分解结果，后续求解同类方程组直接复用），之后每个b只需要简单的点积运算，处理大量b时的总计算量会从O(kn²)降到O(n³ + kn)（k是b的数量），当k很大时，后者的优势极其明显。

补充注意事项 #

如果你的A不是对称矩阵，只需要把预计算步骤改成求解A' \ e1即可，求解器会自动适配转置矩阵的稳定分解，精度和效率不受影响。

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