Menyelesaikan sistem linear

Sistem persamaan linear adalah bentuk masalah yang paling sering dijumpai pada komputasi teknik, dan MATLAB menyediakan satu operator untuk semuanya, yaitu backslash pada pernyataan x = A\b. Topik ini membahas penggunaan operator tersebut untuk sistem persegi yang jumlah persamaannya sama dengan jumlah variabelnya, disertai penyusunan matriks dari masalah rangkaian listrik sederhana. Alasan menghindari penyelesaian lewat inv(A)*b dibahas dari sisi kecepatan dan kestabilan numerik. Perilaku operator yang sama pada sistem over-determined turut dibahas, termasuk solusi kuadrat terkecil yang dihasilkan. Setelah membaca topik ini, sistem linear bentuk apa pun dapat diselesaikan dengan satu idiom baku beserta pemeriksaan mutu hasilnya.

Bentuk masalah dan operator backslash#

Sistem linear ditulis sebagai perkalian matriks,

dengan $A$ menyimpan koefisien, $\mathbf{x}$ vektor variabel yang dicari, dan $\mathbf{b}$ vektor konstanta. Operator backslash A\b menyelesaikan persamaan itu tanpa menampilkan langkah antara apa pun, dan nama resminya adalah mldivide. Operator ini tidak menjalankan satu algoritma saja, melainkan memeriksa struktur $A$ lebih dahulu, kemudian memilih dekomposisi yang paling sesuai. Contohnya, substitusi langsung digunakan untuk matriks segitiga, dekomposisi LU dengan pivoting untuk matriks persegi umum, dan dekomposisi QR untuk matriks persegi panjang. Satu operator yang menutup seluruh kasus tersebut menjadi alasan backslash disebut sebagai idiom baku penyelesaian sistem pada MATLAB.

Pemeriksaan mutu hasil dilakukan dengan mengalikan kembali solusi ke persamaan asalnya. Selisih $A\mathbf{x} - \mathbf{b}$ disebut sisa, dan nilainya yang kecil menandakan solusi memenuhi persamaan hingga batas pembulatan bilangan floating point. Sisa yang besar biasanya menandakan matriks nyaris tunggal, dan penyebabnya dapat diperiksa lebih lanjut dengan rank atau cond dari matriks koefisien.

Contoh: arus pada rangkaian resistor#

Rangkaian listrik sederhana memberi contoh sistem tiga persamaan yang lahir langsung dari hukum fisika. Sebuah sumber 12 V mengalirkan arus $i_1$ melalui resistor 2 Ω, kemudian arus bercabang menjadi $i_2$ melalui resistor 3 Ω dan $i_3$ melalui resistor 6 Ω, dua cabang yang dipasang sejajar. Hukum arus pada simpul percabangan menyatakan arus masuk sama dengan arus keluar. Kedua cabang sejajar tersebut memiliki tegangan yang sama besar, dan hukum tegangan pada loop sumber menyatakan jumlah tegangan sama dengan tegangan sumber,

Ketiga persamaan itu ditata menjadi matriks koefisien $(3 \times 3)$ dan vektor kolom konstanta, lalu diselesaikan dengan satu operator.

% sistem_arus.m: menyelesaikan tiga persamaan arus rangkaian dengan A\b
%
% i1 - i2 - i3 = 0      hukum arus pada simpul percabangan
% i2 - 2*i3    = 0      tegangan kedua cabang sejajar sama besar
% 2*i1 + 3*i2  = 12     hukum tegangan pada loop sumber

A = [ 1  -1  -1;
      0   1  -2;
      2   3   0];
b = [0; 0; 12];

x = A \ b;
fprintf('i1 = %.1f A\n', x(1));
fprintf('i2 = %.1f A\n', x(2));
fprintf('i3 = %.1f A\n', x(3));

disp('sisa A*x - b:')
disp(A * x - b)

i1 = 3.0 A
i2 = 2.0 A
i3 = 1.0 A
sisa A*x - b:
  -2.2204e-16
  -4.4409e-16
            0

Solusinya adalah $i_1 = 3$ A, $i_2 = 2$ A, dan $i_3 = 1$ A, angka bulat yang mudah diperiksa ulang dengan tangan. Sisa $A\mathbf{x} - \mathbf{b}$ berorde $10^{-16}$, yaitu skala pembulatan bilangan floating point, sehingga solusi memenuhi ketiga persamaan sampai presisi mesin. Pemeriksaan sisa seperti ini sebaiknya menjadi kebiasaan penutup setiap penyelesaian sistem karena biayanya satu baris sedangkan informasinya besar, terutama ketika koefisien berasal dari pengukuran yang tidak bulat.

Mengapa bukan inv(A)*b#

Cara klasik pada buku teks, yaitu $\mathbf{x} = A^{-1}\mathbf{b}$, dapat ditulis sebagai inv(A)*b dan memang menghasilkan solusi yang sama secara teori. Praktiknya berbeda pada dua sisi. Pertama, biaya komputasi membentuk invers penuh jauh lebih besar daripada dekomposisi yang hanya memecahkan sistem untuk satu $\mathbf{b}$, sehingga A\b lebih cepat. Kedua, langkah tambahan pada pembentukan invers memperbesar error pembulatan, sehingga A\b juga lebih presisi. Dokumentasi resmi MathWorks untuk fungsi inv menyatakan hal yang sama, yaitu penyelesaian lewat pembagian matriks lebih baik dari sisi waktu maupun akurasi dibandingkan invers eksplisit. Perbedaan presisinya terlihat nyata pada matriks nyaris tunggal seperti matriks Hilbert berorde sepuluh berikut.

H = hilb(10);
b = ones(10, 1);
x_bksl = H \ b;
x_invers = inv(H) * b;
fprintf('kondisi H     : %.3e\n', cond(H));
fprintf('sisa H\\b      : %.3e\n', norm(H * x_bksl - b));
fprintf('sisa inv(H)*b : %.3e\n', norm(H * x_invers - b));

kondisi H     : 1.603e+13
sisa H\b      : 1.769e-10
sisa inv(H)*b : 3.783e-05

Matriks Hilbert memiliki entri $1/(i+j-1)$ yang membuat kolom-kolomnya nyaris sejajar, sehingga kondisinya mencapai $1{,}6 \times 10^{13}$. Pada kondisi setinggi itu, A\b menghasilkan sisa $1{,}8 \times 10^{-10}$ sedangkan inv(A)*b menghasilkan sisa $3{,}8 \times 10^{-5}$, lima orde lebih besar. Kedua nilai sisa inilah ukuran yang tepat untuk dibandingkan karena keduanya deterministik, tidak seperti waktu eksekusi yang berubah antar komputer. Untuk matriks sehat seperti matriks rangkaian sebelumnya, kedua cara memberi hasil yang sama pada tampilan biasa, sehingga alasan memilih backslash tidak terlihat sampai kasus sulit datang.

Sistem over-determined dan kuadrat terkecil#

Pengukuran berkali-kali menghasilkan lebih banyak persamaan daripada variabel, misalnya enam pasang data kalibrasi untuk dua parameter garis lurus. Sistem semacam ini disebut over-determined dan umumnya tidak memiliki solusi persis karena data jarang jatuh tepat pada satu garis. Operator backslash tetap dapat digunakan, dan untuk matriks persegi panjang berpangkat penuh, A\b mengembalikan solusi kuadrat terkecil, yaitu vektor yang memperkecil norma sisa

Contoh berikut mencocokkan garis kalibrasi sensor tekanan melalui enam pengukuran. Matriks rancangan $X$ berukuran $(6 \times 2)$ memuat tekanan pada kolom pertama dan angka satu pada kolom kedua. Pada model garis lurus, kolom pertama dikalikan kemiringan dan kolom kedua dikalikan ofset, sehingga solusinya berisi kedua parameter garis sekaligus.

% regresi_sensor.m: mencocokkan garis kalibrasi sensor lewat kuadrat terkecil
%
% enam pengukuran kalibrasi: keluaran sensor (mV) pada tekanan (bar)
% model garis lurus: keluaran = kemiringan*tekanan + ofset

tekanan   = [0 1 2 3 4 5]';                 % bar
keluaran  = [0.0 2.1 4.2 6.0 8.1 10.2]';    % mV

X = [tekanan, ones(6, 1)];                  % matriks rancangan 6 x 2
fprintf('rank X : %d dari %d baris\n', rank(X), size(X, 1));

sol = X \ keluaran;
fprintf('kemiringan : %.4f mV/bar\n', sol(1));
fprintf('ofset      : %.4f mV\n', sol(2));

sisa = X * sol - keluaran;
fprintf('norm sisa  : %.4f mV\n', norm(sisa));

rank X : 2 dari 6 baris
kemiringan : 2.0229 mV/bar
ofset      : 0.0429 mV
norm sisa  : 0.1757 mV

Nilai rank(X) sama dengan jumlah kolomnya, yaitu dua, sehingga kedua parameter ditentukan secara unik oleh data. Garis hasilnya adalah keluaran $\approx 0.043 + 2.023 \times$ tekanan dengan norma sisa 0.176 mV yang mencerminkan simpangan pengukuran dari garis. Apabila kolom-kolom $X$ sejajar, rank menjadi kurang dari jumlah kolomnya, solusi tidak lagi unik, dan MATLAB maupun Octave mengeluarkan peringatan bahwa matriks kekurangan pangkat. Pemeriksaan rank sebelum menyelesaikan karena itu murah dan langsung menjelaskan status data yang sedang dicocokkan.

Kesalahan umum#

• Menulis b sebagai vektor baris. Backslash mengharapkan ruas kanan sebagai kolom yang tiap barisnya satu persamaan, sehingga b berupa vektor baris ditolak dengan error dimensi. Tulis b = [0; 0; 12], bukan b = [0 0 12].

A = [ 1 -1 -1; 0 1 -2; 2 3 0];
b = [0 0 12];           % vektor baris, bukan kolom
try
  x = A \ b
catch e
  disp(e.message)
end

operator \: nonconformant arguments (op1 is 3x3, op2 is 1x3)

• Menyelesaikan lewat inv(A)*b karena terbiasa dari notasi buku teks. Notasi $A^{-1}\mathbf{b}$ adalah pernyataan teori, sedangkan perhitungan numerik cukup memfaktorkan $A$ sekali untuk satu $\mathbf{b}$. Gunakan A\b sebagai bentuk baku, dan bentuk inv(A) hanya untuk ketika inversnya memang dibutuhkan sebagai objek tersendiri.
• Mengabaikan peringatan matriks nyaris tunggal. Peringatan bahwa matriks tunggal atau kekurangan pangkat berarti solusi yang dikembalikan tidak dapat dipercaya pada digit-digit terakhirnya. Periksa rank(A) dan cond(A) ketika peringatan itu muncul, lalu pertimbangkan mengubah model atau mengumpulkan data yang lebih membedakan variabelnya.

Ringkasan#

• Operator backslash A\b menyelesaikan $A\mathbf{x} = \mathbf{b}$ dengan memilih dekomposisi yang sesuai struktur matriksnya, dan menjadi idiom baku untuk sistem linear.
• Mutu solusi diperiksa dengan sisa $A\mathbf{x} - \mathbf{b}$ yang seharusnya berorde pembulatan mesin.
• inv(A)*b lebih lambat dan kurang presisi daripada A\b, dan perbedaannya menjadi besar pada matriks berkondisi tinggi seperti matriks Hilbert.
• Untuk sistem over-determined berpangkat penuh, A\b menghasilkan solusi kuadrat terkecil sehingga regresi linier dapat dikerjakan tanpa fungsi khusus.

Referensi#

• mldivide, \ — MathWorks — Rujukan resmi operator backslash beserta aturan pemilihan algoritmanya.
• inv — MathWorks — Rujukan invers yang memuat perbandingan resmi antara A\b dan inv(A)*b.
• rank — MathWorks — Rujukan fungsi rank untuk memeriksa keunikan solusi.