0

I'm hoping I'm just missing a simple trick of matrix arithmetic, but the issue I'm having is that all I have access to is an LU solver (Matlab LU* or SuperLU) and I need an LDL decomposition of a symmetric matrix A. So I thought "no problem", since an LU decomposition is unique and an LDL decomposition is unique then D is just the diagonals of U.

But the issue is that it never does A = LU, but rather PA = LU and PA isn't symmetric! So I can't figure out how to determine A = LDL from PA = LU

Is there something easy I can do to fix this? Any help is greatly appreciated.

P.S. Pre-emptive: Yes, I really, really, really, really do need an LDL decomposition. Not there's not another option. Yes, I'm sure. No, you don't need me to lay out the exact problem for you to confirm that I do indeed need an LDL decomposition.

*I know Matlab has its own LDL function, but I'm only using it to prototype and the libraries available to me in C++ (mostly SuperLU) don't seem to have any such function.

cantgetno197
  • 1
  • 1
  • LDL factorization also involves permutations and might give block diagonal D. So I don't quite understand what you want to achieve here. You can use Eigen for extensive C matrix computation libs. Or call DSYTRF from LAPACK directly. – percusse Oct 12 '17 at 12:20
  • I guess the "pre-emptive" was insufficiently pre-emptive. The exact thing I'm trying to accomplish is this exact algorithm: https://dl.acm.org/citation.cfm?id=1916464 I'm writing a solver for electron structure codes. I have no means through which to add or change the libraries included in the software I'm helping to develop. That's out of my hands. – cantgetno197 Oct 12 '17 at 12:31
  • Then you have to implement the Bunch-Kaufman algorithm or something newer yourself. That's a bizarre specification. – percusse Oct 12 '17 at 12:34
  • Ya, I know. I definitely agree. It's a phenomenal pain that can be super frustrating when making quantum transport solvers and especially trying to make use of existing libraries. You really, really, really, need the inverse of a Hermitian indefinite matrix (or at least the diagonals of the inverse). Yes, really, really, really. And the mantra of numerics outside quantum transport is "No one needs an inverse. If you think you need an inverse you must not understand your problem." – cantgetno197 Oct 12 '17 at 12:50
  • Are you then looking for the trace of the inverse? That mantra is always correct since inverse is also `A\eye(n)` – percusse Oct 12 '17 at 13:07
  • No. I mean, this has been a topic of research in computational physics for the last 40 years, if there was a simple way of not doing the inverse, I feel like someone would have figured it out by now. The diagonal of the inverse is the charge density as a function of position (i.e. Ainv(i,i) is the charge at location i). – cantgetno197 Oct 12 '17 at 13:19
  • I'm pretty sure that they won't calculate the whole inverse just for the diagonals. It sounds like an XY problem. – percusse Oct 12 '17 at 13:21
  • No, over the years there have been a number of methods developed to allow you to only calculate the diagonal values of the inverse: SelInv,SINV,RGF, HSC. SelInv is what I'm implementing right now. It needs an LDLT decomposition – cantgetno197 Oct 12 '17 at 13:25
  • If you're really curious I can give you some papers, but the gist is that the matrix is sparse, Hermitian and indefinite and you can get the diagonals without getting the full inverse under some set of conditions (uniform mesh, rectangular system, etc.). I'm trying to make a solver that relaxes some of those conditions but I don't have full control over the libraries available. – cantgetno197 Oct 12 '17 at 13:34
  • "It sounds like an XY problem." Yes, 100% of people on forums like this think this, and I guess, that physicists don't know how math works. Which is why I tried to do the pre-emptive, which never works. – cantgetno197 Oct 12 '17 at 13:58
  • It is not about the preemptive part. What I mean is that probably people much more knowledgeable than me came up with partial solvers. We have similar problems regarding hermitian inverses and there are many modifications to these algorithms especially in sparse linalg domain. You might probably get a better answer from https://scicomp.stackexchange.com/ – percusse Oct 12 '17 at 14:01
  • Well the issue is that I'm doing methods research. There are a small handful of quantum transport solvers out there but they all rely on non-ideal physical assumptions. I'm working on a new type of solver, but it needs to fit into an existing "transport simulation framework" (for funding reasons), and I'm 90% there, I just need to get an LDLT decomposition. – cantgetno197 Oct 12 '17 at 14:09
  • Honestly I didn't understand anything about the specs or your needs. But again LDLt is also gives the permuted result no different than LU. But anyway you can get an almost decomposition with doing another LU on the U matrix of the previous factorization. Or just copy paste an open-source version of the LDL algorithm into your project. I don't see how a framework is an issue to a math library but anyways. Good luck. – percusse Oct 12 '17 at 15:38

1 Answers1

1

Use [L,U] = LU_decomposition(A). Then compute D=Uinv(transpose(L)). Then it holds A = LU = L * D * L'. Notice that D is upper right triangular since inv(transpose(L)) and U are upper right triangular.

Obviously, if A is symmetric then D is symmetric and you have your LDLT decomposition. Hope that helps.

Cheers, Martin

Martin
  • 11
  • 1