分解复数环上矩阵的 LDL 分解

如何解决分解复数环上矩阵的 LDL 分解

我正在使用 Eigen 库进行 LDL 分解。但是，我不是在复数（或实数）环上工作，而是在 typescript primitive string type 环上工作。

我在下面使用 Cholesky-Crout 算法编写了自己的 LDL 分解版本，并对其进行了测试以证明它有效：

template<int M,int N,typename T>
tuple<Matrix<T,M,N>,Matrix<T,N>> ldlt(Matrix<T,N> m)
{
    Matrix<T,N> L = Matrix<T,N>::Identity();
    Matrix<T,N> D;
    for(int i = 0; i < M; i++) {
        for(int j = 0; j < N; j++) {
            if (i == j) {
                T the_sum = 0;
                for(int k = 0; k < j; k++) {
                    the_sum += L(j,k) * conj(L(j,k)) * D(k,k);
                }
                D(j,j) = m(j,j) - the_sum;
            } else if(i > j) {
                T the_sum = 0;
                for(int k = 0; k < j; k++) {
                    the_sum += L(i,k);
                }
                L(i,j) = 1/D(j,j) * (m(i,j) - the_sum);
            }
        }
    }
    return {L,D};
}

此外，我的版本也适用于复数。

我已经尝试使用 Eigen 的 LDLT 版本，但是当我将分解中的矩阵（包括置换矩阵）相乘时，我没有取回我开始使用的矩阵。我的代码如下：

//------------------------------------------------------------------------------
#include <iostream>
#include <Eigen/Core>
#include <Eigen/Dense>
#include <Eigen/Cholesky>
//------------------------------------------------------------------------------
using namespace std;
using namespace Eigen;
//------------------------------------------------------------------------------
#define SPLITCPXCHAR 'j'
//------------------------------------------------------------------------------
class SplitComplex
{
    double re;
    double im;
public:
    SplitComplex();
    SplitComplex(const SplitComplex&); // copy constructor
    SplitComplex(double,double);
    SplitComplex(double);
    double real(void) const;
    double imag(void) const;
    double abs(void) const;
};
//------------------------------------------------------------------------------
// arithmetic operators
SplitComplex operator+(SplitComplex,SplitComplex);
SplitComplex operator-(SplitComplex,SplitComplex);
SplitComplex operator*(SplitComplex,SplitComplex);
SplitComplex operator/(SplitComplex,SplitComplex);
//------------------------------------------------------------------------------
SplitComplex::SplitComplex()
{
    re = 0.0;
    im = 0.0;
}
//------------------------------------------------------------------------------
SplitComplex::SplitComplex(const SplitComplex& aComplex)
{
    re = aComplex.real();
    im = aComplex.imag();
}
//------------------------------------------------------------------------------
SplitComplex::SplitComplex(double d)
{
    re = d;
    im = 0.0;
}
//------------------------------------------------------------------------------
SplitComplex::SplitComplex(double areal,double aimag)
{
    re = areal;
    im = aimag;
}
//------------------------------------------------------------------------------
double SplitComplex::real(void) const
{
    return re;
}
//------------------------------------------------------------------------------
double SplitComplex::imag(void) const
{
    return im;
}
//------------------------------------------------------------------------------
double SplitComplex::abs(void) const
{
    return sqrt(std::abs(re*re-im*im));
}
//------------------------------------------------------------------------------
ostream& operator<<(ostream& aStream,SplitComplex aComplex)
{
    aStream << aComplex.real();
    aComplex.imag() < 0.0 ?
        aStream << " - " << SPLITCPXCHAR << "*"<< -aComplex.imag() :
            aStream << " + " << SPLITCPXCHAR << "*" <<  aComplex.imag();
    return aStream;
}
//------------------------------------------------------------------------------
// arithmetic functions - NOT member of SplitComplex
SplitComplex operator+(SplitComplex aComplex1,SplitComplex aComplex2)
{
    double creal = aComplex1.real() + aComplex2.real();
    double cimag = aComplex1.imag() + aComplex2.imag();
    return SplitComplex(creal,cimag);
}
//------------------------------------------------------------------------------
// arithmetic functions - NOT member of SplitComplex
SplitComplex operator-(SplitComplex aComplex1,SplitComplex aComplex2)
{
    double creal = aComplex1.real() - aComplex2.real();
    double cimag = aComplex1.imag() - aComplex2.imag();
    return SplitComplex(creal,cimag);
}
//------------------------------------------------------------------------------
// arithmetic functions - NOT member of SplitComplex
SplitComplex operator*(SplitComplex aComplex1,SplitComplex aComplex2)
{
    double creal = aComplex1.real() * aComplex2.real() + aComplex1.imag() * aComplex2.imag();
    double cimag = aComplex1.real() * aComplex2.imag() + aComplex1.imag() * aComplex2.real();
    return SplitComplex(creal,cimag);
}
//------------------------------------------------------------------------------
// arithmetic functions - NOT member of SplitComplex
SplitComplex operator/(SplitComplex aComplex1,SplitComplex aComplex2)
{
    double creal = aComplex1.real() * aComplex2.real() - aComplex1.imag() * aComplex2.imag();
    double cimag = -aComplex1.real() * aComplex2.imag() + aComplex1.imag() * aComplex2.real();
    creal /= aComplex2.real() * aComplex2.real() - aComplex2.imag() * aComplex2.imag();
    cimag /= aComplex2.real() * aComplex2.real() - aComplex2.imag() * aComplex2.imag();
    return SplitComplex(creal,cimag);
}
//------------------------------------------------------------------------------
SplitComplex& operator+=(SplitComplex& w,SplitComplex z)
{
    w = w + z;
    return w;
}
//------------------------------------------------------------------------------
SplitComplex& operator-=(SplitComplex& w,SplitComplex z)
{
    w = w - z;
    return w;
}
//------------------------------------------------------------------------------
SplitComplex& operator*=(SplitComplex& w,SplitComplex z)
{
    w = w * z;
    return w;
}
//------------------------------------------------------------------------------
SplitComplex& operator/=(SplitComplex& w,SplitComplex z)
{
    w = w / z;
    return w;
}
//------------------------------------------------------------------------------
bool operator==(SplitComplex w,SplitComplex z)
{
    return w.imag() == z.imag() && w.real() == z.real();
}
//------------------------------------------------------------------------------
namespace Eigen {
    template<> struct NumTraits<SplitComplex>
    : NumTraits<complex<double>>
    {
    typedef double Real;
    typedef SplitComplex NonInteger;
    typedef SplitComplex nested;
    enum {
        IsComplex = 1,IsInteger = 0,IsSigned = 1,RequireInitialization = 1,Readcost = 2,Addcost = 2,MulCost = 4
    };
    };
}
//------------------------------------------------------------------------------
// Declare the split-complex unit
const SplitComplex j(0,1);
//------------------------------------------------------------------------------
// Makes a hermitian split-complex matrix out of an arbitrary real matrix
template<int M,int N>
Matrix<SplitComplex,N> hermitianize(Matrix<double,N>& m)
{
    return m * (1 + j)/2 + m.transpose() * (1 - j)/2;
}
//------------------------------------------------------------------------------
template<int M,D};
}
//------------------------------------------------------------------------------
inline SplitComplex conj(SplitComplex x)  { return x.real() - j * x.imag(); }
inline double real(SplitComplex x)  { return x.real(); }
inline double imag(SplitComplex x)  { return x.imag(); }
inline double abs(SplitComplex x)  { return x.abs(); }
inline double abs2(SplitComplex x)  { return x.abs()*x.abs(); }
inline SplitComplex sqrt(SplitComplex z)
{
    double x = z.real();
    double y = z.imag();
    return (sqrt(x + y) + sqrt(x - y))/2 + j * (sqrt(x + y) - sqrt(x - y))/2;
}
//------------------------------------------------------------------------------
int main(void)
{
    Matrix<double,5,5> m = Matrix<double,5>::Random();
    Matrix<SplitComplex,5> h = hermitianize(m);
    cout << "The split-complex hermitian matrix h is:" << endl << h << endl << endl; // should be a hermitian split-complex matrix
    LDLT<Matrix<SplitComplex,5>> bla = h.ldlt();
    Matrix<SplitComplex,5> L = bla.matrixL();
    Matrix<SplitComplex,5> D = bla.vectorD().asDiagonal();
    Transpositions<5,5> P = bla.transpositionsP();
    cout << "The error from Eigen's LDLT is:" << endl << P.inverse() * L * D * L.adjoint() * P - h << endl << endl; // should be 0 matrix
    auto [LL,DD] = ldlt(h);
    cout << "The error from my LDLT is:" << endl << LL * DD * LL.adjoint() - h << endl << endl; // should be 0 matrix
    return 0;
}
//------------------------------------------------------------------------------