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.
$\begingroup$
$\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
$\begingroup$
$\endgroup$
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.
-
$\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
