共轭梯度法C++实现收敛至错误值问题排查请求

我用C++实现了一个简单的共轭梯度法函数，逻辑上并不复杂，但它无法正确求解线性方程组——虽然会收敛，但结果完全错误。以下是完整代码（默认ApplyLHS和innerProduct函数工作正常）：

共轭梯度法实现代码 #

template <typename TData>
void DoConjGrad(
    std::vector<TData> &in,
    std::vector<TData> &out,
    std::vector<TData> &A
)
{
    TData tol = 10e-8;

    size_t N = in.size();

    std::vector<TData> r(N);
    std::vector<TData> p(N);
    std::vector<TData> w(N);
    size_t k = 0;

    TData rTr;
    TData rTr_prev;
    TData pTw;

    TData alpha;
    TData beta;

    // initial guess for x0
    std::fill(out.begin(), out.end(), 0);

    // r0 = b - Ax0
    ApplyLHS(out, r, A);
    std::transform(in.cbegin(), in.cend(), r.cbegin(), r.begin(), 
        [](const TData &bElem, const TData &rElem){ return bElem - rElem; } );

    // p0 := r0
    std::copy(r.cbegin(), r.cend(), p.begin());

    // calculate (rTr)0
    rTr = innerProduct(r, r);

    for (int i = 0; i < 100; ++i)
    {
        ApplyLHS(p, w, A);
        pTw = innerProduct(p, w);

        alpha = rTr / pTw;

        std::transform(out.cbegin(), out.cend(), p.cbegin(), out.begin(),
            [alpha](const TData &xElem, const TData &pElem)
            {
                return xElem + alpha*pElem;
            });

        std::transform(r.cbegin(), r.cend(), w.cbegin(), r.begin(),
            [alpha](const TData &rElem, const TData &wElem)
            {
                return rElem - alpha*wElem;
            });

        if (rTr < tol)
            break;

        rTr_prev = rTr;
        rTr = innerProduct(r, r);

        beta = rTr / rTr_prev;

        std::transform(r.cbegin(), r.cend(), p.cbegin(), p.begin(),
            [beta](const TData &rElem, const TData &pElem)
            {
                return rElem + beta*pElem;
            });
    }
}

Main函数代码 #

int main()
{
    srand(100);
    size_t N = 2;

    std::vector<double> A(N * N);
    std::vector<double> M(N * N);
    std::vector<double> x_true(N);
    std::vector<double> x_calc(N);
    std::vector<double> b(N);
    
    for (int i = 0; i < N; ++i)
        for (int j = 0; j <= i; ++j)
        {
            double r = 2*((double)rand())/((double)RAND_MAX) - 1;
            M[i*N + j] = r;
            M[j*N + i] = r;
        }
    
    // get positive semi-definite matrix = M.T @ M
    for (int i = 0; i < N; ++i)
        for (int j = 0; j < N; ++j)
            for (int k = 0; k < N; ++k)
                A[i*N + j] = M[k*N + i] * M[k*N + j];

    // generate random x
    for (int i = 0; i < N; ++i)
        x_true[i] = 2*((double)rand())/((double)RAND_MAX) - 1;

    // calculate b
    ApplyLHS(x_true, b, A);

    // solve for x calc
    DoConjGrad<double>(b, x_calc, A);
}

辅助函数代码 #

template <typename TData>
void printMat(
    size_t nr,
    size_t nc,
    std::vector<TData> &mat
)
{
    for (int i = 0; i < nr; ++i)
    {
        for (int j = 0; j < nc; ++j)
        {
            std::cout.precision(5);
            std::cout.width(10);
            std::cout << mat[i*nc + j] << "\t";
        }
        std::cout << "\n";
    }
}

template <typename TData>
void ApplyLHS(
    std::vector<TData> &in,
    std::vector<TData> &out,
    std::vector<TData> &A
)
{
    size_t N = in.size();
    std::fill(out.begin(), out.end(), 0);

    for (int i = 0; i < N; ++i)
    for (int j = 0; j < N; ++j)
        out[i] += in[j] * A[i*N + j];
}

template <typename TData>
TData innerProduct(
    std::vector<TData> &a,
    std::vector<TData> &b
)
{
    TData result = 0;
    for (int i = 0; i < a.size(); ++i)
    {
        result += a[i] * b[i];
    }
    return result;
}

问题根源及修正方案 #

1. 矩阵A初始化错误（最核心问题）

在生成正定矩阵A = M^T @ M时，代码直接赋值而非累加，导致A矩阵仅保留了最后一次k循环的计算结果，完全不是正确的正定矩阵。

修正：
先将A初始化为全0，再执行累加操作：

// 初始化A为全0
std::fill(A.begin(), A.end(), 0.0);
// get positive semi-definite matrix = M.T @ M
for (int i = 0; i < N; ++i)
    for (int j = 0; j < N; ++j)
        for (int k = 0; k < N; ++k)
            A[i*N + j] += M[k*N + i] * M[k*N + j];

2. 收敛判断时机错误

在更新残差r后，代码使用更新前的rTr（残差平方）做收敛判断，导致迭代可能提前终止或无法正确判断收敛状态。

修正：
先计算更新后的残差平方rTr，再进行收敛判断：

// 更新x和r之后
rTr_prev = rTr;
rTr = innerProduct(r, r);

// 移到此处判断收敛
if (rTr < tol)
    break;

beta = rTr / rTr_prev;
// 后续更新搜索方向p

3. 可选调试优化

可以在迭代循环中添加残差和计算结果的输出，方便跟踪收敛过程：

// 在循环内添加
std::cout << "Iteration " << i << ", residual: " << sqrt(rTr) << "\n";
std::cout << "Current x: ";
for (auto val : out) std::cout << val << " ";
std::cout << "\n";

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