Pemecahan sistem linear tercepat untuk matriks persegi kecil (10x10)

9

Saya sangat tertarik untuk mengoptimalkan pemecahan sistem linear untuk matriks kecil (10x10), kadang-kadang disebut matriks kecil . Apakah ada solusi siap untuk ini? Matriks dapat diasumsikan nonsingular.

Pemecah ini harus dijalankan lebih dari 1.000.000 kali dalam mikrodetik pada CPU Intel. Saya berbicara dengan tingkat pengoptimalan yang digunakan dalam permainan komputer. Tidak masalah jika saya mengkodekannya dalam perakitan dan khusus arsitektur, atau mempelajari pengurangan pengorbanan presisi atau keandalan dan menggunakan hack floating point (saya menggunakan flag kompilasi -fast-matematika, tidak ada masalah). Penyelesaiannya bahkan bisa gagal sekitar 20% dari waktu!

ParsialPivLu Eigen adalah yang tercepat di benchmark saya saat ini, mengungguli LAPACK ketika dioptimalkan dengan -O3 dan kompiler yang baik. Tapi sekarang saya pada titik kerajinan tangan pemecah linear kustom. Saran apa pun akan sangat dihargai. Saya akan membuat solusi saya open source dan saya akan mengetahui wawasan kunci dalam publikasi, dll.

Terkait: Kecepatan menyelesaikan sistem linear dengan matriks blok diagonal Apa metode tercepat untuk membalikkan jutaan matriks? https://stackoverflow.com/q/50909385/1489510

rfabbri
sumber
7
Ini terlihat seperti tujuan peregangan. Mari kita asumsikan kita menggunakan Skylake-X Xeon Platinum 8180 tercepat dengan throughput puncak teoritis dari 4 TFLOP presisi tunggal, dan bahwa satu sistem 10x10 membutuhkan sekitar 700 (kira-kira 2n ** 3/3) operasi titik-mengambang yang harus diselesaikan. Kemudian batch 1M sistem tersebut secara teoritis dapat diselesaikan dalam 175 mikrodetik. Itu adalah angka kecepatan cahaya yang tidak bisa melebihi. Bisakah Anda membagikan kinerja apa yang saat ini Anda capai dengan kode tercepat yang ada? BTW, apakah data presisi tunggal atau presisi ganda?
njuffa
@ njuffa ya saya bertujuan untuk mencapai hampir 1 ms tetapi mikro adalah cerita lain. Untuk mikro saya mempertimbangkan mengeksploitasi struktur invers tambahan dalam batch dengan mendeteksi matriks yang serupa, yang sering terjadi. Perf saat ini berada pada kisaran 10-500 ms tergantung pada prosesor. Presisi ganda atau bahkan kompleks ganda. Presisi tunggal memang lebih lambat.
rfabbri
@njuffa Saya dapat mengurangi atau meningkatkan presisi untuk kecepatan
rfabbri
2
Sepertinya presisi / ketepatan bukanlah prioritas Anda. Untuk tujuan Anda, mungkin metode berulang yang terpotong pada sejumlah kecil evaluasi berguna? Apalagi jika Anda memiliki tebakan awal yang masuk akal.
Spencer Bryngelson
1
Apakah Anda berputar? Bisakah Anda melakukan faktorisasi QR bukan eliminasi Gaussian. Apakah Anda interleave sistem Anda sehingga Anda dapat menggunakan instruksi SIMD dan melakukan beberapa sistem sekaligus? Apakah Anda menulis program garis lurus tanpa loop dan tanpa pengalamatan tidak langsung? Akurasi apa yang Anda inginkan dan bagaimana saya akan mengkondisikan sistem Anda? Apakah mereka memiliki struktur yang dapat dieksploitasi.
Carl Christian

Jawaban:

7

Menggunakan tipe matriks Eigen di mana jumlah baris dan kolom dikodekan ke dalam tipe pada waktu kompilasi memberi Anda keunggulan atas LAPACK, di mana ukuran matriks hanya diketahui saat runtime. Informasi tambahan ini memungkinkan kompiler untuk melakukan loop penuh atau sebagian membuka gulungan, menghilangkan banyak instruksi cabang. Jika Anda ingin menggunakan pustaka yang ada daripada menulis kernel Anda sendiri, memiliki tipe data di mana ukuran matriks dapat dimasukkan sebagai parameter template C ++ mungkin akan sangat penting. Satu-satunya perpustakaan lain yang saya tahu yang melakukan ini adalah blaze , sehingga mungkin layak dilakukan benchmark terhadap Eigen.

Jika Anda memutuskan untuk menggulung implementasi Anda sendiri, Anda mungkin menemukan apa yang dilakukan PETSc untuk format CSR bloknya menjadi contoh yang berguna, meskipun PETSc sendiri mungkin tidak akan menjadi alat yang tepat untuk apa yang ada dalam pikiran Anda. Daripada menulis loop, mereka menulis setiap operasi tunggal untuk perkalian matriks-vektor kecil secara eksplisit (lihat file ini di repositori mereka). Ini menjamin bahwa tidak ada instruksi cabang seperti yang Anda dapatkan dengan satu loop. Versi kode dengan instruksi AVX adalah contoh yang baik tentang bagaimana sebenarnya menggunakan ekstensi vektor. Misalnya, fungsi ini menggunakan__m256dtipe data untuk beroperasi secara bersamaan pada empat ganda pada saat yang sama. Anda bisa mendapatkan peningkatan kinerja yang cukup besar dengan menuliskan secara eksplisit semua operasi menggunakan ekstensi vektor, hanya untuk faktorisasi LU alih-alih penggandaan matriks-vektor. Daripada benar-benar menulis kode C dengan tangan, Anda akan lebih baik menggunakan skrip untuk menghasilkannya. Mungkin juga menyenangkan untuk melihat apakah ada perbedaan kinerja yang cukup besar ketika Anda memesan ulang beberapa operasi untuk lebih memanfaatkan perpipaan instruksi.

Anda mungkin juga mendapatkan beberapa jarak tempuh dari alat STOKE , yang secara acak akan mengeksplorasi ruang kemungkinan transformasi program untuk menemukan versi yang lebih cepat.

Daniel Shapero
sumber
tx. Saya sudah menggunakan Eigen seperti Map <const Matrix <complex, 10, 10>> AA (A) berhasil. akan memeriksa ke dalam hal-hal lain.
rfabbri
Eigen juga memiliki AVX dan bahkan header.h kompleks untuk itu. Mengapa PETSc untuk ini? Sulit untuk bersaing dengan Eigen dalam hal ini. Saya mengkhususkan Eigen bahkan lebih untuk masalah saya dan dengan strategi pivot perkiraan yang bukannya mengambil max di atas kolom, menukar pivot segera ketika menemukan yang lain yaitu 3 kali lipat lebih besar.
rfabbri
1
@ rfabbri Saya tidak menyarankan Anda menggunakan PETSc untuk ini, hanya saja apa yang mereka lakukan dalam contoh khusus itu bisa menjadi pelajaran. Saya telah mengedit jawaban untuk membuatnya lebih jelas.
Daniel Shapero
4

Gagasan lain dapat menggunakan pendekatan generatif (program menulis program). Penulis program (meta) yang mengeluarkan urutan instruksi C / C ++ untuk melakukan unpivoted ** LU pada sistem 10x10 .. pada dasarnya mengambil sarang loop k / i / j dan meratakannya menjadi O (1000) atau lebih baris dari aritmatika skalar. Kemudian masukkan program yang dihasilkan ke mana pun yang mengoptimalkan kompiler. Apa yang saya pikir agak menarik di sini, adalah menghapus loop memperlihatkan setiap ketergantungan data dan redeksan subekspresi, dan memberikan kompiler kesempatan maksimum untuk menyusun ulang instruksi sehingga mereka memetakan dengan baik ke perangkat keras yang sebenarnya (misalnya jumlah unit eksekusi, bahaya / kios, sehingga di).

Jika Anda mengetahui semua matriks (atau bahkan hanya beberapa dari mereka), Anda dapat meningkatkan throughput dengan memanggil intrinsik / fungsi SIMD (SSE / AVX) alih-alih kode skalar. Di sini Anda akan mengeksploitasi paralelisme yang memalukan di seluruh instance, alih-alih mengejar paralelisme apa pun dalam satu instance. Sebagai contoh, Anda dapat melakukan 4 presisi ganda LU secara simultan menggunakan intrinsik AVX256, dengan mengemas 4 matriks "melintasi" register dan melakukan operasi yang sama ** pada semuanya.

** Karenanya fokus pada LU yang tidak diproteksi. Berputar merusak pendekatan ini dalam dua cara. Pertama, ini memperkenalkan cabang karena pemilihan pivot, artinya dependensi data Anda tidak begitu dikenal. Kedua, ini berarti bahwa "slot" SIMD yang berbeda harus melakukan hal-hal yang berbeda, karena instance A mungkin berputar berbeda dari instance B. Jadi, jika Anda mengejar semua ini, saya sarankan secara statis memutar matriks Anda sebelum perhitungan (permutasi entri terbesar dari setiap kolom ke diagonal).

rchilton1980
sumber
karena matriksnya sangat kecil, mungkin pivoting dapat dihilangkan jika pra-diskalakan. Bahkan tidak melakukan pra-putar matriks. Yang kita butuhkan hanyalah entri berada dalam 2-3 urutan besarnya satu sama lain.
rfabbri
2

Pertanyaan Anda mengarah pada dua pertimbangan berbeda.

Pertama, Anda harus memilih algoritma yang tepat. Oleh karena itu, pertanyaan apakah matriks memiliki struktur, harus dipertimbangkan. Misalnya, ketika matriks simetris, dekomposisi Cholesky lebih efisien daripada LU. Ketika Anda hanya membutuhkan keakuratan terbatas, metode berulang bisa lebih cepat.

Kedua, Anda perlu mengimplementasikan algoritma secara efisien. Untuk melakukannya, Anda perlu mengetahui hambatan dari algoritma Anda. Apakah implementasi Anda terikat oleh kecepatan transfer memori atau oleh kecepatan perhitungan. Karena Anda hanya mempertimbangkan matriks, matriks Anda harus masuk ke dalam cache CPU sepenuhnya. Dengan demikian, Anda harus menggunakan unit SIMD (SSE, AVX, dll.) Dan inti prosesor Anda, untuk melakukan sebanyak mungkin komputasi per siklus.10×10

Secara keseluruhan, jawaban atas pertanyaan Anda sangat tergantung pada perangkat keras dan matriks yang Anda pertimbangkan. Mungkin tidak ada jawaban yang pasti dan Anda harus mencoba beberapa hal untuk menemukan metode yang optimal.

H. Rittich
sumber
Sejauh ini Eigen sudah sangat mengoptimalkan, menggunakan SEE, AVX, dll dan saya mencoba metode iteratif dalam tes pendahuluan dan mereka tidak membantu. Saya mencoba Intel MKL tetapi tidak lebih baik dari Eigen dengan flag GCC yang dioptimalkan. Saat ini saya mencoba untuk membuat sesuatu yang lebih baik dan lebih sederhana daripada Eigen dan untuk melakukan tes yang lebih rinci dengan metode berulang.
rfabbri
1

Saya akan mencoba inversi blockwise.

https://en.wikipedia.org/wiki/Invertible_matrix#Blockwise_inversion

Eigen menggunakan rutin yang dioptimalkan untuk menghitung kebalikan dari matriks 4x4, yang mungkin merupakan yang terbaik yang akan Anda dapatkan. Coba gunakan itu sebanyak mungkin.

http://www.eigen.tuxfamily.org/dox/Inverse__SSE_8h_source.html

Kiri atas: 8x8. Kanan atas: 8x2. Kiri bawah: 2x8. Kanan bawah: 2x2. Balikkan 8x8 menggunakan kode inversi 4x4 yang dioptimalkan. Sisanya adalah produk matriks.

EDIT: Menggunakan blok 6x6, 6x4, 4x6, dan 4x4 terbukti sedikit lebih cepat daripada yang saya jelaskan di atas.

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;
}

Berikut adalah hasil dari satu run mark bench menggunakan satu juta Eigen::Matrix<double,10,10>::Random()matriks dan Eigen::Matrix<double,10,1>::Random()vektor. Dalam semua pengujian saya, kebalikan saya selalu lebih cepat. Pemecahan rutin saya melibatkan menghitung invers dan kemudian mengalikannya dengan vektor. Terkadang lebih cepat dari Eigen, terkadang tidak. Metode menandai bangku saya mungkin cacat (tidak menonaktifkan turbo boost, dll). Juga, fungsi acak Eigen mungkin tidak mewakili data nyata.

  • Pembalikan poros parsial eigen: 3036 milidetik
  • Pembalikan saya dengan blok atas 8x8: 1638 milidetik
  • Pembalikan saya dengan blok atas 6x6: 1234 milidetik
  • Pivot parsial Eigen teratasi: 1791 milidetik
  • Solusi saya dengan blok atas 8x8: 1739 milidetik
  • Solusi saya dengan blok atas 6x6: 1286 milidetik

Saya sangat tertarik untuk melihat apakah ada yang bisa mengoptimalkan ini lebih lanjut, karena saya memiliki aplikasi elemen hingga yang membalikkan matriks 10x10 trilyun (dan ya, saya memang membutuhkan koefisien individu dari invers sehingga secara langsung menyelesaikan sistem linear tidak selalu merupakan pilihan) .

Charlie S
sumber