1
$\begingroup$

I wanted to generate LU decomposition of large size dense matrices ($N>10^7$), the LU decomposition I'm currently using is based on Adaptive Cross Approximation and is taking very long time to execute for larger $N$. Can anybody suggest a few LU decomposition techniques that can be very well parallelised (using OpenMP) and take a shorter period of time?

Note:

  1. I write the code in C++ and make use of Xeon Processor(128 threads) and Eigen library.
  2. The entries in the matrix are filled through a kernel function of form $e^{-(x_1-x_2)^2}$.
  3. Storage of matrix is not a problem, I'm working on Xeon processor and has enough memory and moreover, I'm not storing the full matrix and whenever I need to find an entry in matrix, I use the kernel function and generate type-double number for that cell.
Improve this question
$\endgroup$
2
  • $\begingroup$ Do you mean your matrix is ~3k*3k or ~10m*10m when you say dense matrix of size N>10^7 ? $\endgroup$
    – DrHansGruber
    Commented Jan 24, 2017 at 15:51
  • $\begingroup$ @DrHansGruber If I take N = 6553600 the matrix is NxN matrix where if u want a cell entry u need to use the kernel, hence not storing the matrix anywhere. $\endgroup$
    – user7440094
    Commented Jan 24, 2017 at 18:31

2 Answers 2

Reset to default
1
$\begingroup$

You can't. Unless you have some special knowledge, your $L$ and $U$ factors are going to be dense, and will have $N^2 = {\cal O}(10^{14})$ entries -- more than you can store on any reasonable machine (about a million GB). Furthermore, it takes ${\cal O}(N^3)={\cal O}(10^{21})$ operations to compute these factors, i.e., $3\cdot 10^4$ CPU years -- longer than you want to wait.

The problem you have can simply not be solved via LU decomposition. You need to come up with other strategies for doing what you want to do -- e.g., solving linear systems with Krylov space methods, computing explicit Green's functions, or using knowledge about sparsity.

Improve this answer
$\endgroup$
4
  • $\begingroup$ I have already run my code for N=6553600 and it works, I repeat again I didn't store the NxN matrix rather get the entries from kernel function.The LU decomposition of such N is possible and gets stored with my processor. $\endgroup$
    – user7440094
    Commented Jan 24, 2017 at 18:29
  • $\begingroup$ @user7440094 -- I do not believe you. If your matrix is $N\times N$, then so are the LU factors, and you need to store $N^2$ floating point values for the LU factors -- which in your case would be somewhere around $4\cdot 10^{13}$, or 40,000 GB. I don't believe you that your methods works because I don't believe that you have that much memory. $\endgroup$
    – Wolfgang Bangerth
    Commented Jan 24, 2017 at 22:40
  • $\begingroup$ Wolfgang, he is not really talking about a full-blown LU decomposition, but an approximation to it using some sort of sparsification technique, probably of hierarchical structure. Such methods can reduce the memory to $N \log N$ if properly implemented (of course, problem dependent). Certainly, it is going to give an approximate LU decomposition, not an exact one. $\endgroup$
    – Anton Menshov ♦
    Commented Jan 25, 2017 at 22:14
  • 1
    $\begingroup$ @AntonMenshov, Hm, I see. The question does not make that clear. $\endgroup$
    – Wolfgang Bangerth
    Commented Jan 26, 2017 at 17:24
1
$\begingroup$

I'm not exactly sure what type of framework for LU are you using, as one can apply ACA to various different setups. And to make it explicit, the approach you are working on right now, and I am going to propose does not give you an LU-decomposition. It is an approximation of an LU-decomposition in some form.

Since you are using a non-singular kernel, I guess you can try to compress the whole matrix. Usually, those methods are applied to singular kernels and then you have a hierarchical pattern, when you subdivide you matrix into blocks in a multilevel fashion and most of them are compressible. Those blocks would be computed using some fast technique, like ACA. It's worth mentioning, that ACA is less than ideal, has troubles with the controlled accuracy and etc. Then, you can factorize the matrix. The approach I described is pretty much Hierarchical matrix ($\mathcal H$-matrix) approach, developed by Dr. W. Hackbusch http://www.hmatrix.org/

To cut along story short, that approach allows for $\mathcal O(k_\text{max}^3N\log^2N)$ LU decomposition and $\mathcal O(k_\text{max}^3N\log N)$ for back substitution and pretty well parallelizable via OpenMP. Notice, here $k_\text{max}$ is the maximum rank, that for ACA would mean the number of skeletons. Of course, that assumes that you have a nice balanced partitioning. Details and limitations are available in numerous papers on $\mathcal H$-matrix papers.

Take a look at this approach, and their available $\mathcal H$lib library, that is available to the public. Other suggestion might include several other frameworks, like HSS, that I know about (but much less familiar), and more from computational electromagnetics:

  1. J. Shaeffer, "Million plus unknown MOM LU factorization on a PC," 2015 International Conference on Electromagnetics in Advanced Applications (ICEAA), Turin, 2015, pp. 62-65. doi: 10.1109/ICEAA.2015.7297075

  2. S. Kapur and D. E. Long, "N-body problems: IES3: Efficient electrostatic and electromagnetic simulation," in IEEE Computational Science and Engineering, vol. 5, no. 4, pp. 60-67, Oct.-Dec. 1998. doi: 10.1109/MCSE.1998.7102081

Improve this answer
$\endgroup$
3
  • $\begingroup$ yes I'm having an approximate of LU decomposition for HODLR( hierarchically off-diagonal low-rank) type matrices , thanks for the resources, I'll read and try to apply the method mentioned. $\endgroup$
    – user7440094
    Commented Jan 25, 2017 at 4:46
  • $\begingroup$ The approach you mention have higher order as compared to ACA which is O(r^2(m+n)) where r is low rank of matrix.So, I think ACA will still be better. $\endgroup$
    – user7440094
    Commented Feb 13, 2017 at 6:46
  • $\begingroup$ I doubt that is the case. As far as I know, you would require that complexity just to approximate the matrix. Computing its approx. LU would require another step that would probably have $k^3$ dependence. $\endgroup$
    – Anton Menshov ♦
    Commented Feb 14, 2017 at 2:25

Your Answer

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

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