1
$\begingroup$

I have implemented a schur complement domain decomposition method for solving large scale problems in Matlab. It works well but I think there should be some alternative approaches to implement some parts of the algorithm more efficiently to circumvent high memory usage. Specifically, I need to evaluate $M=BQ^{-1}B^T$ where $Q$ is a $p\times p$ real-valued, symmetric, sparse and positive-definite matrix and $B$ is a $q\times p$ extremely sparse matrix only containing $\{-1,+1\}$ nonzero entries ($q$ is usually much greater than $p$). My approach is $M=B*(Q\backslash B^T)$; however its evaluation becomes very memory- and time-consuming for not very large values of $p$. Any suggestion to get around this problem is warmly welcome.

$\endgroup$
  • $\begingroup$ Upon closer scrutiny, do you perhaps mean $p$ is much greater than $q$? (such that $\mathbf M$ is smaller than $\mathbf Q$)? If $\mathbf M$ is in fact larger than $\mathbf Q$, I think you have no choice but to leave it in that "outer product" form and apply it as a cascade of 3 operations. $\endgroup$ – rchilton1980 Mar 2 '18 at 3:48
  • $\begingroup$ Or perhaps you meant to say $\mathbf Q$ was size $q \times q$ and $\mathbf B$ is size $p \times q$ with $q>>p$? $\endgroup$ – rchilton1980 Mar 2 '18 at 3:50
  • $\begingroup$ No, $q>>p$ is correct according to the given dimensions in the question. $\endgroup$ – Mat123 Mar 2 '18 at 7:24
  • $\begingroup$ What do $\mathbf Q$ and $\mathbf B$ represent? My initial guess was that $\mathbf Q$ was a domain solve, and $\mathbf B$ was a surface restriction/projection, but the indicated sizes suggest otherwise. Perhaps $\mathbf Q$ is instead a coarse space solve, and $\mathbf B$ represents the effect of updating every domain with a coarse space correction? $\endgroup$ – rchilton1980 Mar 5 '18 at 14:18
2
$\begingroup$

This is always going to be a somewhat costly calculation because the inverse of any "interesting" sparse matrix is generally dense and therefore so is $\mathbf M$. That said, there are smarter techniques than the obvious one you mention, that are asymptotically no slower than the factorization of $\mathbf Q$ itself, as long as $\mathbf B$ has no more rows than the root separator of $\mathbf Q$. This is likely the case, because in most DDM's the matrix $\mathbf Q$ comes from a volume discretization of the underlying PDE over a single domain (volume) and the operator $\mathbf M$ represents a Dirichlet-to-Neumann map or Robin-to-Robin map upon the boundary of a domain, such that the $\mathbf M$ has "surface cardinality" just like the root separator of $\mathbf Q$.

I wrote a GPL-licensed sparse direct package called Myramath to perform exactly this calculation (see http://www.myracore.com, it's also in my profile). In addition to the usual factor()/solve() API, the solver also possesses a method called .schur() that can form a matrix like $\mathbf M$. I'd also recommend another algorithm called .partialsolve(), which also has interesting applications for substructuring/DDM. See the tutorial material for details, it starts here. The schur() and partialsolve() methods are addressed here.

$\endgroup$
  • $\begingroup$ Thanks! Unfortunately I do not think it is relevant for your problem. The schur() method is useful for "condensation-like" operations, where high dimensional PDE's are projected upon low dimensional boundaries. From our discussion about the relative sizes/shapes of $\mathbf Q$ and $\mathbf B$, you apparently have the opposite case. $\endgroup$ – rchilton1980 Mar 2 '18 at 12:22
  • $\begingroup$ You say that "there are smarter techniques", but you do not mention or describe them at all. As it is, this is more or less a link-only answer. $\endgroup$ – Federico Poloni Mar 3 '18 at 11:03
  • 1
    $\begingroup$ The smarter technique: given symmetric positive definite $\mathbf Q$, compute the Cholesky factor $\mathbf {Q=LL^T}$, then the schur complement $\mathbf {M = B*Q^{-1}B'=[L^{-1}B^T]^T[L^{-1}B^T]}$. The key insight is that $\mathbf{L^{-1}B^{T}}$ is sparse, and can be calculated using sparsity-aware/multifrontal methods just like the original factorization itself. With some effort, this can be extended to other flavors of matrices too (non-positive, asymmetric, etc). Details (and source code) are in the links. $\endgroup$ – rchilton1980 May 1 '18 at 13:22
  • $\begingroup$ That said, this doesn't appear to address the OP's question anyway, and they seem to have lost interest in the problem. $\endgroup$ – rchilton1980 May 1 '18 at 13:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.