Ah good hint. Apparently a slightly more general version is called the *-Sylvester Equation
http://gauss.uc3m.es/web/personal_web/fteran/papers/equation_lama_rev.pdf and quite some work has been done, though I haven't time to dig through the references. The big difference in your case is that L and M are lower triangular, and L appears in both terms multiplying the unknown, so there's a lot more structure than usual (which leads to the fact that forward substitution works). ~Chris On Tue, Jun 3, 2014 at 1:17 PM, Thomas Covert <[email protected]> wrote: > Thanks for that, Chris. I also worked out a similar derivation - though I > wasn't able to prove to myself that the equation B = L*M' + M*L' has a > unique solution for M. > > This feels similar to the Sylvester Equation, but its not quite the same... > > > On Mon, Jun 2, 2014 at 10:13 PM, Chris Foster <[email protected]> wrote: >> >> On Tue, Jun 3, 2014 at 4:43 AM, Thomas Covert <[email protected]> >> wrote: >> > I was hoping to find some neat linear algebra trick that would let me >> > compute a DualNumber cholesky factorization without having to resort to >> > non-LAPACK code, but I haven't found it yet. That is, I figured that I >> > could compute the cholesky factorization of the real part of the matrix >> > and >> > then separately compute a matrix which represented the DualNumber part. >> > However, I've not yet found a math text describing such a beast. >> >> I'm not sure exactly where you'd look this up, but here's a go at >> setting the problem up. >> >> Suppose you're trying to factorise >> >> (A + e*B) = (L + e*M) * (L + e*M)' >> >> where e is the dual imaginary unit. That is, (L + e*M) is the desired >> Cholesky >> factor of (A + e*B). Expanding this out using e*e = 0 leads to >> >> A = L*L' >> B = L*M' + M*L' >> >> So you can compute the real part L of the Cholesky decomposition exactly >> as >> usual. Given that you now have L, you want the lower triangular matrix M. >> Because L and M are lower triangular that's actually quite easy: matrix >> elements of M can be computed from the top left to bottom right in a >> process >> that basically looks like forward substitution. >> >> You probably want this second step to use LAPACK though, since it's still >> O(N^3)... Ideally the system of equations >> >> B = L*M' + M*L' >> >> actually has a name and there's a standard way to solve it. I don't know >> it >> though and a quick search failed to enlighten me. Perhaps someone else >> can >> give us a hint? >> >> ~Chris > >
