# Efficient evaluation of $BQ^{-1}B^T$ (Domain Decomposition Implementation)

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.

• 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. – rchilton1980 Mar 2 '18 at 3:48
• Or perhaps you meant to say $\mathbf Q$ was size $q \times q$ and $\mathbf B$ is size $p \times q$ with $q>>p$? – rchilton1980 Mar 2 '18 at 3:50
• No, $q>>p$ is correct according to the given dimensions in the question. – Mat123 Mar 2 '18 at 7:24
• 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? – rchilton1980 Mar 5 '18 at 14:18

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.
• 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. – rchilton1980 Mar 2 '18 at 12:22
• 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. – rchilton1980 May 1 '18 at 13:22