问题

我希望解决一般的线性方程组A*x = b.m-by-m矩阵是稀疏的,实数,正方形,条件稍差,非对称,但它是奇异的(秩(A)== m-1),因为x只知道一个加性常数.

我可以通过在三个矢量(指定其非零项创建矩阵A i,j和v),使得A(i(k),j(k)) = v(k):

A = sparse( i, j, v, m, m );

原始方程

我可以解决这个原始方程如下:

x = A \ b;

如果我想要一个独特的解决方案,我可以在计算非唯一解决方案后施加一个约束(比如,x(4)== 3.14159):

x = x - x(4) + 3.14159;

修正方程

我可以C通过添加一个额外的唯一性约束来创建一个新的全秩矩阵,如下所示:

% Add the constraint x(4) == 3.14159
extraRow = zeros(1,m);
extraRow(4) = 1.0;
C = [A; extraRow];    % Add to the matrix A
d = [b; 3.14159];     % Add to the RHS vector, …

任何人都可以给我一个简单的例子,说明如何使用CHOLMOD将元素添加到三元组矩阵中.

我尝试过这样的事情:

cholmod_triplet *A;
int k;

void add_A_entry(int r, int c, double x)
{
    ((int*)A->i)[k] = r;
    ((int*)A->j)[k] = c;
    ((double*)A->x)[k] = x;
    k++;
}

int main()
{
    k = 0;
    cholmod_common com;
    cholmod_start(&com);

    A = cholmod_allocate_triplet(202, 202, 202*202, -1, CHOLMOD_REAL, &com);
    add_A_entry(2, 2, 1.);
    add_A_entry(4, 1, 2.);
    add_A_entry(2, 10, -1.);

    cholmod_print_triplet(A, "A", &com);

    cholmod_finish(&com);
    return 0;
}

但是,这不会向矩阵添加任何元素.我只是得到输出:

CHOLMOD triplet: A:  202-by-202, nz 0, lower.  OK

当然,我试图通过搜索和CHOLMOD文档找到解决方案,但我没有找到任何帮助.

我目前正在做一个我们需要解决的项目

|Ax - b|^2.

在这种情况下，A是一个非常稀疏的矩阵A'A，每行最多有 5 个非零元素。

我们正在处理图像，其中的维度A'A是NxNN 是像素数。在这种情况下N = 76800。我们计划去RGB然后维度将是3Nx3N。

在 matlab 中(A'A)\(A'b)，使用双精度求解大约需要 0.15 秒。

我现在已经对Eigens稀疏求解器进行了一些试验。我试过了：

SimplicialLLT
SimplicialLDLT
SparseQR
ConjugateGradient

和一些不同的排序。迄今为止最好的是

SimplicialLDLT

这需要0.35 - 0.5使用AMDOrdering.

例如，当我使用ConjugateGradient它时，大约需要6 s，0用作初始化。

解决问题的代码是：

    A_tot.makeCompressed();
     // Create solver
    Eigen::SimplicialLDLT<Eigen::SparseMatrix<float>, Eigen::Lower, Eigen::AMDOrdering<int> > solver;
    // Eigen::ConjugateGradient<Eigen::SparseMatrix<float>, Eigen::Lower> cg;
    solver.analyzePattern(A_tot);
    t1 = omp_get_wtime();
    solver.compute(A_tot);
    if (solver.info() != Eigen::Success) …

在SparseSuiteQR中，我可以找到的所有示例都使用 stdin 或读取的文件来创建稀疏矩阵。有人可以提供一个如何直接在 C++ 中创建一个简单的示例吗？

更好的是，在 CHOLMOD 文档中，提到了 matlab 中可用的稀疏2函数，其行为与稀疏函数相同。这可以在C++中使用吗？