增广矩阵的行减少 - 3D 样条计算

Question

我正在尝试通过松弛三次样条实现插值，可在本文第 5 章（第 9 页）中找到： https://www.math.ucla.edu/~baker/149.1.02w/handouts/dd_splines.pdf

到目前为止，我有以下内容：

auto GetControlPoints = [](const std::vector<Vector3d>& S) {
    int n = S.size();
    float var = n - 1.0f;
    MatrixXd M(n - 1, n - 1);
    VectorXd C[3] = {
        VectorXd(n - 1),
        VectorXd(n - 1),
        VectorXd(n - 1)
    };

    for (int i = 0; i < n - 1; ++i) {
        auto r = RowVectorXd(n - 1);

        for (int j = 0; j < n - 1; ++j) {
            if (j == i)
                r[j] = var;
            else if (j == i - 1 || j == i + 1)
                r[j] = 1.f;
            else
                r[j] = 0.f;
        }

        M.row(i) = r;

        if (i == 0) {
            for (int j = 0; j < 3; ++j) {
                C[j] << (n + 1) * S[1][j] - S[0][j];
            }
        }
        else if (i == n - 1) {
            for (int j = 0; j < 3; ++j) {
                C[j] << (n + 1) * S[n - 1][j] - S[n][j];
            }
        }
        else {
            for (int j = 0; j < 3; ++j) {
                C[j] << (n + 1) * S[i][j];
            }
        }
    }

    MatrixXd augMC[3] = {
        MatrixXd(n - 1, n),
        MatrixXd(n - 1, n),
        MatrixXd(n - 1, n)
    };

    for (int i = 0; i < 3; ++i) {
        augMC[i].block(0, 0, n - 1, n - 1) = M;
        augMC[i].block(n - 1, n - 1, n - 1, 1) = C[i].transpose();
    }
};

我已经到了使用 M 和 C 制作增广矩阵的地步，但我不知道如何对其进行行缩减。有什么想法吗？

Answer 1

您可以使用 inplace-variant of PartialPivLU——但看起来您实际上想要求解系统 M*B = C，您应该为其分解 M（因为它是对称的，您可以使用 LLt 或 LDLt 分解），然后使用该分解的 solve 方法。

要设置 M，您还应该使用 diagonal 方法（未测试）：

MatrixXd M(n - 1, n - 1);
M.setZero();
M.diagonal().setConstant(n - 1.0);
M.diagonal<1>().setOnes();
M.diagonal<-1>().setOnes();

LLT<MatrixXd> lltOfM(M);
for (int i = 0; i < 3; ++i) { B[i] = lltOfM.solve(C[i]); }

对于大 n 这是次优的，因为它没有利用 M 的三对角结构。您可以为此尝试稀疏模块，但实际上应该有一个直接算法（尽管 Eigen 没有明确具有三对角矩阵类型）。

对于 C，您可能还可以使用 MatrixX3d（我真的不明白您是如何填充 C 向量的——我认为您当前的代码应该在 运行-time，除非n==2，或者你禁用了断言）。