Performa pada Kode Ilmiah

Kode ilmiah yang menghasilkan jawaban benar belum tentu memenuhi tenggat komputasinya. Topik ini membahas tiga pengungkit performa yang paling sering menentukan pada kode ilmiah Julia, yaitu pemilihan solver yang sesuai sifat masalah, pemanfaatan sparsity pada sistem besar, dan pengendalian alokasi memori pada fungsi dinamika. Prinsip stabilitas tipe dari modul Performa dan Stabilitas Tipe menjadi latar bagi ketiga pembahasan tersebut.

Memilih solver sesuai sifat masalah#

Solver eksplisit seperti Tsit5 menghitung langkah berikutnya langsung dari evaluasi fungsi dinamika, sehingga setiap langkahnya murah. Keuntungan itu hilang pada persamaan kaku (stiff), yaitu masalah yang memuat komponen dengan skala waktu yang berbeda jauh, misalnya peluruhan dengan laju 1000 per satuan waktu yang diganggu gaya luar yang berubah perlahan. Pada masalah seperti itu, batas stabilitas memaksa solver eksplisit memperkecil langkah jauh di bawah kebutuhan akurasinya, sehingga jumlah langkah membengkak.

Solver implisit bekerja dengan cara berbeda: setiap langkah menyelesaikan persamaan nonlinear yang melibatkan Jacobian $\partial f/\partial u$. Satu langkah menjadi lebih mahal, tetapi panjang langkahnya tidak dibatasi stabilitas. Keluarga Rodas termasuk metode Rosenbrock yang menjadi pilihan umum untuk masalah kaku berukuran kecil hingga menengah, dan contoh berikut menggunakan Rodas5P sebagai anggotanya.

using OrdinaryDiffEq

# Masalah kaku: laju peluruhan 1000 kali lebih cepat daripada variasi gaya luar
f_stiff!(du, u, p, t) = (du[1] = -1000.0 * (u[1] - sin(t)))

masalah = ODEProblem(f_stiff!, [1.0], (0.0, 10.0))

solusi_eksplisit = solve(masalah, Tsit5())
solusi_implicit = solve(masalah, Rodas5P())

println("Tsit5   : ", length(solusi_eksplisit.t), " langkah")
println("Rodas5P : ", length(solusi_implicit.t), " langkah")
println("u akhir Tsit5   : ", round(solusi_eksplisit.u[end][1], digits = 5))
println("u akhir Rodas5P : ", round(solusi_implicit.u[end][1], digits = 5))

Tsit5   : 2863 langkah
Rodas5P : 31 langkah
u akhir Tsit5   : -0.54306
u akhir Rodas5P : -0.54318

Perbandingan 2.863 langkah dengan 31 langkah memperlihatkan gejala khas masalah kaku. Kedua solver menghasilkan lintasan yang sama-sama dapat diterima pada toleransi bawaannya, tetapi biaya langkah eksplisit yang murah tidak menutupi ledakan jumlah langkahnya. Aturan praktisnya, mulailah dengan Tsit5 dan pindah ke solver implisit apabila jumlah langkah membengkak tanpa kebutuhan akurasi yang setara. Nilai u akhir yang berbeda pada angka keempat desimal juga menunjukkan bahwa kedua solver mengikuti toleransi bawaannya masing-masing, bukan menghitung lintasan yang identik.

Kode pada contoh sekaligus memperlihatkan manfaat antarmuka bersama. Objek masalah yang sama diberikan kepada dua solver yang berbeda tanpa perubahan apa pun pada fungsi dinamikanya, dan OrdinaryDiffEq menyediakan ratusan metode lain untuk kebutuhan di antara keduanya. Percobaan solver semacam ini murah karena hanya argumen solve yang berubah, sehingga pemilihan solver dapat dilakukan berdasarkan pengukuran, bukan dugaan.

Memanfaatkan sparsity#

Sistem ilmiah berukuran besar sering menghasilkan Jacobian yang mayoritas elemennya nol, misalnya karena setiap variabel hanya berinteraksi dengan tetangganya. Menyimpan matriks seperti itu secara padat berarti memboroskan memori dan operasi untuk elemen kosong. Pustaka standar SparseArrays menyimpan hanya elemen non-nol dalam format CSC (compressed sparse column), sehingga ukuran memori dan biaya perkaliannya mengikuti jumlah elemen non-nol, bukan jumlah total sel. Konstruksi spdiagm pada contoh menyusun matriks langsung dari diagonal-diagonalnya, sedangkan fungsi nnz menghitung banyaknya elemen non-nol yang tersimpan.

using SparseArrays, LinearAlgebra

n = 10_000
A = spdiagm(-1 => fill(-1.0, n - 1), 0 => fill(4.0, n), 1 => fill(-1.0, n - 1))

println("Dimensi  : ", size(A))
println("nnz      : ", nnz(A))
println("Kepadatan: ", nnz(A) / (n^2))

b = A * ones(n)
println("b[1] = ", b[1], ", b[2] = ", b[2], ", b[n] = ", b[n])

Dimensi  : (10000, 10000)
nnz      : 29998
Kepadatan: 0.00029998
b[1] = 3.0, b[2] = 2.0, b[n] = 3.0

Matriks tridiagonal berukuran sepuluh ribu kali sepuluh ribu hanya menyimpan 29.998 elemen non-nol, sedangkan versi padatnya memuat 100 juta entri atau sekitar 800 MB untuk elemen Float64. Perkalian matriks-vektor pada bentuk sparse hanya menjumlahkan pasangan yang benar-benar ada, sehingga hasil pada baris pertama, kedua, dan terakhir tetap dihitung dari diagonal 4.0 serta tetangga -1.0 yang diharapkan.

Dalam konteks suite SciML, pola sparsity itu diteruskan ke solver melalui argumen jac_prototype pada ODEProblem. Argumen tersebut menerima matriks sparse sebagai cetakan struktur Jacobian, sehingga faktorisasi yang dilakukan solver implisit memanfaatkan pola nolnya. Untuk sistem besar, langkah ini sering menentukan apakah simulasi dapat diselesaikan sama sekali atau justru terjebak pada faktorisasi matriks yang hampir seluruhnya kosong.

Mengendalikan alokasi pada fungsi dinamika#

Fungsi dinamika adalah kode yang paling sering dieksekusi pada simulasi ODE karena dipanggil untuk setiap evaluasi pada setiap langkah, dan jumlah pemanggilannya mudah mencapai jutaan. Setiap alokasi array baru di dalam fungsi tersebut berulang sebanyak itu pula. Gaya in-place f(du, u, p, t) menghindarinya dengan menulis hasil langsung ke du yang disiapkan solver, sedangkan gaya out-of-place f(u, p, t) membuat array baru pada setiap pemanggilan.

using OrdinaryDiffEq

# In-place: menulis hasil ke du yang sudah tersedia
function dinamika!(du, u, p, t)
    du[1] = u[2]
    du[2] = -p[1] * u[1] - p[2] * u[2]
end

# Out-of-place: mengembalikan array baru pada setiap pemanggilan
dinamika(u, p, t) = [u[2], -p[1] * u[1] - p[2] * u[2]]

function ukur(u, p, du)
    a1 = @allocated dinamika(u, p, 0.0)
    a2 = @allocated dinamika!(du, u, p, 0.0)
    c1 = @allocations dinamika(u, p, 0.0)
    c2 = @allocations dinamika!(du, u, p, 0.0)
    a1, a2, c1, c2
end

p = (1.0, 0.25)
u = [1.0, 0.0]
du = similar(u)

ukur(u, p, du)  # pemanasan pertama agar kode telah terkompilasi
a1, a2, c1, c2 = ukur(u, p, du)
println("Out-of-place: ", a1, " byte, ", c1, " alokasi")
println("In-place    : ", a2, " byte, ", c2, " alokasi")

Out-of-place: 80 byte, 2 alokasi
In-place    : 0 byte, 0 alokasi

Pengukuran dilakukan dengan @allocated, yaitu makro yang menghitung byte yang dialokasikan selama evaluasi ekspresi, disertai @allocations yang menghitung banyaknya kejadian alokasi. Keduanya ditempatkan di dalam fungsi dan setelah pemanasan, karena kompilasi pertama serta cakupan global menambah angka yang tidak mewakili kerja rutinnya. Gaya out-of-place tetap mengalokasikan pada setiap panggilan, sedangkan versi in-place berjalan tanpa alokasi sama sekali. Selisih yang tampak kecil pada satu pemanggilan menjadi besar ketika dikalikan ribuan langkah dan banyak evaluasi per langkah.

Pola ini terkait erat dengan stabilitas tipe. Fungsi dinamika yang tidak stabil tipenya memaksa Julia mengalokasikan nilai yang seharusnya berada di register, dan gejalanya sama persis dengan pengukuran di atas. Argumen dan hasil yang bertipe konkret serta mutasi pada du tanpa membuat array baru menjaga bagian simulasi yang paling sering dieksekusi tetap bebas alokasi.

Ringkasan#

• Persamaan kaku membatasi panjang langkah solver eksplisit melalui stabilitas, sehingga solver implisit dari keluarga Rodas menyelesaikannya dengan jauh lebih sedikit langkah.
• SparseArrays menyimpan hanya elemen non-nol, dan argumen jac_prototype meneruskan pola sparsity Jacobian kepada solver implisit.
• Fungsi dinamika in-place menulis ke du tanpa alokasi, dan @allocated setelah pemanasan memverifikasinya.
• Stabilitas tipe dan pola mutasi tanpa array baru berjalan bersamaan untuk menjaga bagian simulasi yang paling sering dieksekusi tetap bebas alokasi.

Referensi#

• DifferentialEquations.jl Documentation — SciML — Dokumentasi solver, termasuk pedoman pemilihan solver untuk masalah kaku dan tidak kaku.
• SparseArrays — The Julia Language — Pustaka standar matriks jarang, dari konstruksi spdiagm hingga fungsi nnz.