Row Reduction of augmented matrix - 3D Spline calculations

Question

I'm trying to implement Interpolation by relaxed cubic splines, which can be found in the 5th chapter of this article (page 9): https://www.math.ucla.edu/~baker/149.1.02w/handouts/dd_splines.pdf

So far, I have the following:

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();
    }
};

I got to the point where I made an augmented Matrix using M and C, but I have no idea on how to row reduce it. Any thoughts?

Answer 1

You could use the inplace-variant of PartialPivLU -- but it looks like you actually want to solve the system M*B = C for which you should just decompose M (as it is symmetric you can use an LLt or LDLt decomposition) and then use the solve method of that decomposition.

To setup M you should also use the diagonal method (untested):

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]); }

For large n this is sub-optimal, since it does not exploit the tridiagonal structure of M. You could try out the sparse module for this, but there should actually be a direct algorithm (Eigen does not explicitly have a tridiagonal matrix type, though).

For C you probably could also use a MatrixX3d (I don't really understand how you fill your C vectors -- I think your current code should assert at run-time, unless n==2, or you disabled assertions).