1
$\begingroup$

I would like to solve a dense linear system the form in python $$ L\left(\boldsymbol{x}\right):=\left[\gamma^+\left[\boldsymbol{A}+\frac{1}{2}\boldsymbol{B}^{-1}\right] +\gamma^-\left[\boldsymbol{A}-\frac{1}{2}\boldsymbol{B}^{-1}\right]\right]\cdot\boldsymbol{x}=\boldsymbol{b} $$ I thought it is a good idea not to explicitly compute the inverse of $\boldsymbol{B}$. Therefore, I implemented the Operator $L$ to solve $L\left(\boldsymbol{x}\right)=\boldsymbol{b}$ using Krylov methods like cg or gmres. I am using the scipy.sparse.linalg.LinearOperator class for the operator see the docs here. The product $\boldsymbol{B}^{-1}\cdot\boldsymbol{x}$ computed by another Krylov iteration or an L-U decomposition depending on system size.

However, for larger problems I would like to improve the rate of convergence of the outer iteration. I neither have sparse matrices nor an explicit representation of my matrix. Therefore, as far as I know the classical preconditioners for Krylov methods like ilu or Jacobi's method are not applicable.

Are there other methods which can be used? And are there python libraries for these methods?

Improve this question
$\endgroup$

1 Answer 1

Reset to default
1
$\begingroup$

I would try to circumvent the problem entirely by the substitution $x=By$. Then your equation reduces to $$ \left[ (\gamma^+ + \gamma^-) AB + \frac{1}{2} (\gamma^+-\gamma^-) I \right] y = b.$$ If $AB$ is sparse, then I would first attempt a direct solve, and failing that, I would apply GMRES with a sparse preconditioner. If successful, then you can recover $x$ from $y$ using your current strategy.

Any extra information about the properties of the matrices $A$ and $B$ will be vital to the selection of the optimal solution strategy. It would be great to know the physical problem which generated the matrices, their dimensions, sparsity pattern, and any special mathematical properties inherited from the underlying problem.

EDIT: On a desktop I would want to avoid large dense matrices. On a larger parallel machine I would give ScaLAPACK a try as the LU factorization algorithm is built on top of the matrix-matrix multiplication operation. This is one of the few kernels which are compatible with the hardware of today. It has high arithmetic intensity, so the processors can run at a high fraction of their peak flop rate.

That being said, I would also investigate if my dense matrix admits a good sparse approximation. I would treat my matrix as a vector and sort all components by their absolute value. This would immediately reveal if the overwhelming majority of the entries are tiny relative to the others. I would construct a preconditioner for my dense matrix by dropping all entries below a suitable threshold.

Improve this answer
$\endgroup$
6
  • $\begingroup$ Thank you for reply. The equation above is simplified. The linear system arises from a boundary element discretization of an interface problem (Laplace equation). In literature, it says that the problem is elliptic and some matrices are even positive definite. I am now trying the reformulate the problem in block-matrix form to avoid the inverse. However, as you may matrices from BEM-methods are allways dense. Therefore, sparsity patterns should be not help, or? $\endgroup$
    – sebastian_g
    Jun 1, 2016 at 9:18
  • $\begingroup$ @sebastian_g i have added the thoughts that came to mind when I read your comment. $\endgroup$
    – Carl Christian
    Jun 1, 2016 at 18:13
  • $\begingroup$ Is the second approach of introducing a drop tolerance to get a sparse approximation of system matrix the same what the incomplete LU-factorization does? I think I can program this sparse approximation myself and than use scipy.sparse.linalg.splu to get a LU-factorization. scipy uses the SuperLU library. I don't know how efficient it is. But I am not worried about processor loads at this stage. $\endgroup$
    – sebastian_g
    Jun 3, 2016 at 7:07
  • $\begingroup$ ILU can apply a variety of dropping rules. In your case, I would start by dropping all entries on a row which a small relative to the diagonal entry. In truth, I do not know that this will work for you application. While you matrix is dense, i.e. everything is connected to everything, I am gambling that the strength of connection decays with distance and that a good sparse approximation is possible. If the diagonal entry is zero or no good, then base the comparison on the largest element in the row instead. $\endgroup$
    – Carl Christian
    Jun 3, 2016 at 7:29
  • $\begingroup$ Meanwhile I implemented the dropping criterion you suggested and works fine. I can now solve my test case up to machine in reasonable time using gmres. I observed that restart is important in my problem. $\endgroup$
    – sebastian_g
    Jun 7, 2016 at 10:55

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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