计算科学 - 计算稀疏矩阵的零空间 - 吾爱随笔录

我正在计算稀疏矩形的零空间 $m$ X $n$ 矩阵 $A$ ，在哪里 $m$ << $n$ . 我通过计算 QR 分解来做到这一点 $A^T$ 并提取 $n-m$ 结果的最右边的列 $Q$ . 象征性地：

A^{T} = Q R, Q = [Q_{1} Q_{2}], Q_{2} = N u l l s p a c e (A)

$A^T=QR, Q=[Q_1 Q_2] , Q_2 = Nullspace(A)$

使用这种技术，我得到了我想要的结果，但我希望它更快。

这是稀疏模式 $A^T$ 由 Mathematica 生成（2,337 x 36,090；42,900 个非零）：

以及由此产生的零空间（36,090 x 33,753；8,167,054 个非零）

稀疏模式源于有限元问题中两个子域接触的约束。左边的系数是正的，而右边的系数是负的。这是由于符号约定。我使用有限元网格中的节点连接矩阵对每个子域的节点索引进行排序，并将其输入 METIS（这使我的共轭梯度求解速度更快！）。

我主要将 Eigen 和 Intel MKL 用于我的矩阵代数，但也将 SuiteSparse 用于他们更快的 QR 分解例程。尽管对于这个问题并不重要，但我将我的 Eigen 和 SuiteSparse 之间的接口过程包括在内，以防万一它确实重要。这是我计算零空间（内核）的方法 $A$ ：

using namespace Eigen;
typedef SparseMatrix<double,ColMajor,int> Sparse_t;

cholmod_common* cc, Common;
cc = &Common;
cholmod_l_start(cc);

Sparse_t AT; //assume A^T is already defined
//unfortunately necessary copy with SuiteSparse_long indices -- minimal impact on performance
SparseMatrix<double, ColMajor, SuiteSparse_long> tmp(AT);

const int rows = AT.rows();
const int cols = AT.cols();
const int ordering = SPQR_ORDERING_BEST;
const double tol = SPQR_NO_TOL;
const SuiteSparse_long econ = rows;
const int getCTX = 1; 

cholmod_sparse AT_cm = viewAsCholmod(tmp);
cholmod_sparse* I_cm = cholmod_l_speye(rows, rows, CHOLMOD_REAL, cc);
cholmod_sparse* Q_cm = NULL;

//perform QR factorization and request (Q^T * I)^T = Q
SuiteSparseQR<double>(ordering, tol, econ, getCTX, &AT_cm, I_cm, NULL, &Q_cm, NULL, NULL, NULL, NULL, NULL, NULL, cc);
Sparse_t kernel = viewAsEigen<Sparse_t::Scalar, Sparse_t::Options, Sparse_t::Index>(*Q_cm).rightCols(rows - cols);
cholmod_l_free_sparse(&I_cm, &cc);
cholmod_l_free_sparse(&Q_cm, &cc);

根据我对该网站和其他来源的研究，QR 分解似乎是计算稀疏矩阵的零空间的最佳方法。此外，SuiteSparse 的 SPQR 似乎是稀疏矩阵最快和最通用的例程。我已经在 Windows 上编译了所有内容，并将 SuiteSparse 与 Intel MKL（可能与此问题无关）和 Intel TBB 链接起来，并且在使用 SPQR 时确实实现了有限的多线程。

我认为我的代码针对我在六核 CPU 笔记本电脑上选择的平台 (Windows) 进行了优化。我没有太多运气让 SPQR 的 GPU 加速功能工作，但这可能是我的弱 GPU。但是，我真正的问题是是否有更好的方法来解决这个问题。我缺少线性代数技巧吗？也许我可以执行一些矩阵分解 $A^T$ 在执行 QR 之前？

谢谢！（我想知道多久有人删除寒暄）

更新：我使用umfpack先计算LU分解，然后计算QR分解。

P_{r} R A P_{c} = L U, U^{T} = Q R, Q = [Q_{1} Q_{2}], N u l l s p a c e (A) = P_{c} Q_{2}

$P_rRAP_c=LU , U^T=QR , Q = [Q_1Q_2] , Nullspace(A) = P_cQ_2$ 我没有发布额外的代码，因为它是一个令人困惑的 c 风格的混乱！零空间是在很短的时间内计算出来的（包括 LU 分解），但由此产生的稀疏模式非常奇怪，这似乎会减慢后续操作的速度。我的数值结果似乎不受影响。如果

A

$A$ 非常稀疏，我的整体有限元计算的改进可以忽略不计。作为

A

$A$ 变得不那么稀疏（更多约束），性能有非常显着的提高。也许不同的节点排序策略可以解决这个问题。

$U^T 的零空间

与 QR 分解相比，LU 分解是瞬时的。有没有办法直接从分解中计算零空间，形式为 $P_rRAP_c=LU$ ?