#c #matrix #linear-algebra #eigen #eigen3
#c #матрица #линейная алгебра #eigen #eigen3
Вопрос:
Я использую собственную библиотеку для декомпозиции ЛПНП. Однако вместо того, чтобы работать над кольцом комплексных чисел (или действительных чисел) Я работаю над кольцом расщепленных комплексных чисел.
Ниже я запрограммировал свою собственную версию разложения ЛПНП, используя алгоритм Холески-Краута, и протестировал его, чтобы показать, что он работает:
template<int M, int N, typename T>
tuple<Matrix<T,M,N>,Matrix<T,M,N>> ldlt(Matrix<T,M,N> m)
{
Matrix<T,M,N> L = Matrix<T,M,N>::Identity();
Matrix<T,M,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) * conj(L(j,k)) * D(k,k);
}
L(i,j) = 1/D(j,j) * (m(i,j) - the_sum);
}
}
}
return {L, D};
}
Кроме того, моя версия работает и над комплексными числами.
Я пробовал использовать собственную версию 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 SplitComplexamp;); // 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 SplitComplexamp; 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));
}
//------------------------------------------------------------------------------
ostreamamp; operator<<(ostreamamp; 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);
}
//------------------------------------------------------------------------------
SplitComplexamp; operator =(SplitComplexamp; w, SplitComplex z)
{
w = w z;
return w;
}
//------------------------------------------------------------------------------
SplitComplexamp; operator-=(SplitComplexamp; w, SplitComplex z)
{
w = w - z;
return w;
}
//------------------------------------------------------------------------------
SplitComplexamp; operator*=(SplitComplexamp; w, SplitComplex z)
{
w = w * z;
return w;
}
//------------------------------------------------------------------------------
SplitComplexamp; operator/=(SplitComplexamp; w, SplitComplex z)
{
w = w / z;
return w;
}
//------------------------------------------------------------------------------
bool operator==(SplitComplex w, SplitComplex z)
{
return w.imag() == z.imag() amp;amp; 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,M,N> hermitianize(Matrix<double,M,N>amp; m)
{
return m * (1 j)/2 m.transpose() * (1 - j)/2;
}
//------------------------------------------------------------------------------
template<int M, int N, typename T>
tuple<Matrix<T,M,N>,Matrix<T,M,N>> ldlt(Matrix<T,M,N> m)
{
Matrix<T,M,N> L = Matrix<T,M,N>::Identity();
Matrix<T,M,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) * conj(L(j,k)) * D(k,k);
}
L(i,j) = 1/D(j,j) * (m(i,j) - the_sum);
}
}
}
return {L, 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,5>::Random();
Matrix<SplitComplex,5,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,5>> bla = h.ldlt();
Matrix<SplitComplex,5,5> L = bla.matrixL();
Matrix<SplitComplex,5,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;
}
//------------------------------------------------------------------------------
Комментарии:
1. «Я не возвращаю матрицу, с которой начал» — это ожидаемо, поскольку математика с плавающей запятой имеет тенденцию накапливать ошибки. Можете ли вы быть более конкретным? В чем заключается ошибка?
2. @Yakk-AdamNevraumont Когда я вычитаю произведение
P.inverse() * L * D * L.adjoint() * Ph, я не получаю матрицу, все элементы которой маленькие. Они не обязательно должны быть точно равны нулю, но они должны быть очень маленькими. В то время как с моей реализацией я получаю приблизительно нулевую матрицу3. Похоже, что существует проблема с продуктами, включающими транспозиции справа. Если вы вычисляете
PtL = P.transpose() * L;, а затемPtL * D * PtL.adjoint()вы находитесь в пределах машинной точности исходной матрицы. Кроме того, это то,bla.reconstructedMatrix()что вы ожидаете (вам нужно предоставитьoperator!=для этого).4. @chtz Если вы напишете ответ, я приму его