Spróbowałbym odwrócić blokowo.
https://en.wikipedia.org/wiki/Invertible_matrix#Blockwise_inversion
Eigen używa zoptymalizowanej procedury do obliczenia odwrotności macierzy 4x4, co jest prawdopodobnie najlepszą możliwą wartością. Spróbuj użyć tego jak najwięcej.
http://www.eigen.tuxfamily.org/dox/Inverse__SSE_8h_source.html
Lewy górny róg: 8x8. U góry po prawej: 8x2. U dołu po lewej: 2x8. Na dole po prawej: 2x2. Odwróć 8x8 za pomocą zoptymalizowanego kodu inwersyjnego 4x4. Reszta to produkty matrycowe.
EDYCJA: Używanie bloków 6x6, 6x4, 4x6 i 4x4 okazało się nieco szybsze niż to, co opisałem powyżej.
using namespace Eigen;
template<typename Scalar, int tl_size, int br_size>
Matrix<Scalar, tl_size + br_size, tl_size + br_size> blockwise_inversion(const Matrix<Scalar, tl_size, tl_size>& A, const Matrix<Scalar, tl_size, br_size>& B, const Matrix<Scalar, br_size, tl_size>& C, const Matrix<Scalar, br_size, br_size>& D)
{
Matrix<Scalar, tl_size + br_size, tl_size + br_size> result;
Matrix<Scalar, tl_size, tl_size> A_inv = A.inverse().eval();
Matrix<Scalar, br_size, br_size> DCAB_inv = (D - C * A_inv * B).inverse();
result.topLeftCorner<tl_size, tl_size>() = A_inv + A_inv * B * DCAB_inv * C * A_inv;
result.topRightCorner<tl_size, br_size>() = -A_inv * B * DCAB_inv;
result.bottomLeftCorner<br_size, tl_size>() = -DCAB_inv * C * A_inv;
result.bottomRightCorner<br_size, br_size>() = DCAB_inv;
return result;
}
template<typename Scalar, int tl_size, int br_size>
Matrix<Scalar, tl_size + br_size, tl_size + br_size> my_inverse(const Matrix<Scalar, tl_size + br_size, tl_size + br_size>& mat)
{
const Matrix<Scalar, tl_size, tl_size>& A = mat.topLeftCorner<tl_size, tl_size>();
const Matrix<Scalar, tl_size, br_size>& B = mat.topRightCorner<tl_size, br_size>();
const Matrix<Scalar, br_size, tl_size>& C = mat.bottomLeftCorner<br_size, tl_size>();
const Matrix<Scalar, br_size, br_size>& D = mat.bottomRightCorner<br_size, br_size>();
return blockwise_inversion<Scalar,tl_size,br_size>(A, B, C, D);
}
template<typename Scalar>
Matrix<Scalar, 10, 10> invert_10_blockwise_8_2(const Matrix<Scalar, 10, 10>& input)
{
Matrix<Scalar, 10, 10> result;
const Matrix<Scalar, 8, 8>& A = input.topLeftCorner<8, 8>();
const Matrix<Scalar, 8, 2>& B = input.topRightCorner<8, 2>();
const Matrix<Scalar, 2, 8>& C = input.bottomLeftCorner<2, 8>();
const Matrix<Scalar, 2, 2>& D = input.bottomRightCorner<2, 2>();
Matrix<Scalar, 8, 8> A_inv = my_inverse<Scalar, 4, 4>(A);
Matrix<Scalar, 2, 2> DCAB_inv = (D - C * A_inv * B).inverse();
result.topLeftCorner<8, 8>() = A_inv + A_inv * B * DCAB_inv * C * A_inv;
result.topRightCorner<8, 2>() = -A_inv * B * DCAB_inv;
result.bottomLeftCorner<2, 8>() = -DCAB_inv * C * A_inv;
result.bottomRightCorner<2, 2>() = DCAB_inv;
return result;
}
template<typename Scalar>
Matrix<Scalar, 10, 10> invert_10_blockwise_6_4(const Matrix<Scalar, 10, 10>& input)
{
Matrix<Scalar, 10, 10> result;
const Matrix<Scalar, 6, 6>& A = input.topLeftCorner<6, 6>();
const Matrix<Scalar, 6, 4>& B = input.topRightCorner<6, 4>();
const Matrix<Scalar, 4, 6>& C = input.bottomLeftCorner<4, 6>();
const Matrix<Scalar, 4, 4>& D = input.bottomRightCorner<4, 4>();
Matrix<Scalar, 6, 6> A_inv = my_inverse<Scalar, 4, 2>(A);
Matrix<Scalar, 4, 4> DCAB_inv = (D - C * A_inv * B).inverse().eval();
result.topLeftCorner<6, 6>() = A_inv + A_inv * B * DCAB_inv * C * A_inv;
result.topRightCorner<6, 4>() = -A_inv * B * DCAB_inv;
result.bottomLeftCorner<4, 6>() = -DCAB_inv * C * A_inv;
result.bottomRightCorner<4, 4>() = DCAB_inv;
return result;
}
Oto wyniki jednego testu porównawczego z użyciem miliona Eigen::Matrix<double,10,10>::Random()macierzy i Eigen::Matrix<double,10,1>::Random()wektorów. We wszystkich moich testach moja odwrotność jest zawsze szybsza. Moja procedura rozwiązywania polega na obliczeniu odwrotności, a następnie pomnożeniu jej przez wektor. Czasami jest szybszy niż Eigen, czasem nie. Moja metoda oznaczania na ławce może być wadliwa (nie wyłączała turbo boost itp.). Ponadto losowe funkcje Eigen mogą nie reprezentować rzeczywistych danych.
- Odwrotność częściowa własna odwrotna: 3036 milisekund
- Moja odwrotność z górnym blokiem 8x8: 1638 milisekund
- Moja odwrotność z górnym blokiem 6x6: 1234 milisekundy
- Częściowe rozwiązanie osi własnych w rozwiązaniu: 1791 milisekund
- Moje rozwiązanie z górnym blokiem 8x8: 1739 milisekund
- Moje rozwiązanie z górnym blokiem 6x6: 1286 milisekund
Bardzo mnie interesuje, czy ktoś może to dalej zoptymalizować, ponieważ mam aplikację elementów skończonych, która odwraca matryce gazillionów 10x10 (i tak, potrzebuję indywidualnych współczynników odwrotności, więc bezpośrednie rozwiązanie układu liniowego nie zawsze jest opcją) .