7 NumPy ile Lineer Cebir
Bilimsel hesaplamadaki problemlerin büyük çoğunluğu sonunda bir lineer cebir problemine dönüşür. Bir diferansiyel denklemi sayısal olarak çözmek, bir veri kümesine eğri uydurmak ya da bir optimizasyon problemini çözmek çoğu zaman bir \(Ax = b\) sistemini çözmeye ya da bir matrisin özdeğerlerini bulmaya varır. Kâğıt üzerinde \(3 \times 3\) bir sistemi Gauss yöntemiyle birkaç dakikada çözeriz; bilgisayar \(2000 \times 2000\) bir sistemi, bu bölümde ölçeceğimiz gibi, saniyenin onda biri kadar bir sürede çözer.
NumPy’ın numpy.linalg modülü bu hesapları onlarca yıldır geliştirilen ve sınanan LAPACK kütüphanesine devreder. Bu bölümde Lineer Cebir notlarından tanıdığımız işlemleri (matris çarpımı, ters matris, determinant, rank, özdeğerler) NumPy ile yapacağız. Tanımları ve ispatları yeniden anlatmayacağız; gerektiğinde Lineer Cebir notlarına bağlantı vereceğiz. Dizilerin nasıl kurulduğunu NumPy Dizileri, eleman eleman işlemleri ve broadcasting’i Vektörizasyon ve Broadcasting bölümünde gördük.
Asıl dikkatimiz, kayan noktalı sayılarla hesap yapmanın getirdiği farklarda olacak. Tekil bir matrisin determinantı \(0\) yerine \(6{,}7 \cdot 10^{-18}\) çıkabilir; solve hiçbir uyarı vermeden bileşenleri \(10^{15}\) mertebesinde olan anlamsız bir çözüm döndürebilir; ters matrisle çözmek hem daha yavaş hem de daha az doğrudur. Bu olguları açıklayan kavram koşul sayısıdır. Bölümü en küçük kareler problemi ve tekil değer ayrışımıyla kapatacağız.
7.1 Matris Çarpımı ve Temel İşlemler
NumPy’da bir matris, iki eksenli bir dizidir: birinci eksen satırları, ikinci eksen sütunları sayar. Matematikteki \(AB\) çarpımının Python’daki karşılığı @ işlecidir.
import numpy as np
A = np.array([[1, 2],
[3, 4]])
B = np.array([[0, 1],
[1, 0]])
print(A @ B) # matris çarpımı
print(B @ A) # sıra değişince sonuç değişir
print(A * B) # eleman eleman çarpım: matris çarpımı DEĞİLÇıktı:
[[2 1]
[4 3]]
[[3 4]
[1 2]]
[[0 2]
[3 0]]\(B\) matrisi \(A\)’yı sağdan çarpınca \(A\)’nın sütunlarının, soldan çarpınca satırlarının yerini değiştirir; bu yüzden \(AB \ne BA\)’dır. Son satırdaki A * B ise matris çarpımı değildir: karşılıklı elemanları çarpar.
Tanım 7.1 (Matris Çarpımı İşleci) \(A\) dizisinin şekli \((m, n)\), \(B\) dizisinin şekli \((n, p)\) olsun. A @ B ifadesi, elemanları
\[(AB)_{ij} = \sum_{k=1}^{n} a_{ik} b_{kj}\]
olan \((m, p)\) şekilli diziyi verir. \(A\)’nın sütun sayısı \(B\)’nin satır sayısına eşit değilse işlem ValueError hatası verir. np.matmul(A, B) aynı işlemin fonksiyon biçimidir.
Yani @, Lineer Cebir notlarındaki matris çarpımının ta kendisidir. Aşağıdaki şekil bu boyut kuralını gösteriyor: çarpımın her elemanı, \(A\)’nın bir satırı ile \(B\)’nin bir sütununun iç çarpımıdır.
A @ B için boyut kuralı: (3, 4) şekilli A ile (4, 2) şekilli B'nin çarpımı (3, 2) şekillidir. C'nin i. satır, j. sütundaki elemanı, A'nın i. satırı ile B'nin j. sütununun iç çarpımıdır. Bu yüzden A'nın sütun sayısı B'nin satır sayısına eşit olmak zorundadır.A * B, aynı şekilli iki diziyi eleman eleman çarpar. Şekiller farklıysa broadcasting kuralları devreye girer ve hiçbir hata almadan beklenmedik bir sonuç elde edilebilir: örneğin \((3, 3)\) şekilli A ile uzunluğu \(3\) olan v için A * v hata vermez ama \(Av\) değildir. Matematikte \(AB\) yazdığımız her yerde kodda A @ B yazılır. Eski kodlarda görülen np.dot(A, B) iki eksenli dizilerde aynı sonucu verir; np.matrix sınıfı ise artık önerilmez ve yeni kodda kullanılmaz.
Tek eksenli bir dizi ne satır ne de sütun vektörüdür; yalnızca bir sayı listesidir. @ işleci onu bağlama göre yorumlar: A @ v çarpımında v bir sütun vektörü, w @ A çarpımında w bir satır vektörü gibi davranır. Aynı uzunluktaki iki vektörün u @ v çarpımı ise standart iç çarpımdır. Transpoz A.T ile alınır.
import numpy as np
A = np.array([[2, -1, 0],
[1, 3, 4]]) # 2 x 3 matris
v = np.array([1, 0, -2]) # uzunluğu 3 olan vektör
w = np.array([1, 1]) # uzunluğu 2 olan vektör
print(A @ v) # A v: v sütun gibi davranır
print(w @ A) # w^T A: w satır gibi davranır
print(v @ v) # iç çarpım: 1 + 0 + 4
print(A.T.shape) # transpoz
print(A @ A.T) # 2 x 2 ve simetrik
print(np.trace(A @ A.T)) # köşegen toplamı (iz)Çıktı:
[ 2 -7]
[3 2 4]
5
(3, 2)
[[ 5 -1]
[-1 26]]
31\(AA^T\) her zaman simetriktir, çünkü \((AA^T)^T = (A^T)^T A^T = AA^T\)’dir. np.trace köşegen elemanların toplamını, yani izi verir. \(n \times n\) birim matris np.eye(n) ile, bir kare matrisin \(k\). kuvveti np.linalg.matrix_power(A, k) ile kurulur.
Örnek 7.1 (Fibonacci Sayıları ve Matris Kuvveti) Fibonacci sayıları \(F_0 = 0\), \(F_1 = 1\) ve \(F_{n+1} = F_n + F_{n-1}\) ile tanımlanır. \(Q = \begin{pmatrix} 1 & 1 \\ 1 & 0 \end{pmatrix}\) matrisinin kuvvetlerini kullanarak \(F_{10}\) ve \(F_{100}\) sayılarını hesaplayınız.
Çözüm
Önce \(n \ge 1\) için
\[Q^n = \begin{pmatrix} F_{n+1} & F_n \\ F_n & F_{n-1} \end{pmatrix}\]
olduğunu tümevarımla görelim. \(n = 1\) için \(F_2 = 1\), \(F_1 = 1\), \(F_0 = 0\) olduğundan eşitlik doğrudur. Eşitlik \(n\) için doğruysa
\[Q^{n+1} = Q^n Q = \begin{pmatrix} F_{n+1} + F_n & F_{n+1} \\ F_n + F_{n-1} & F_n \end{pmatrix} = \begin{pmatrix} F_{n+2} & F_{n+1} \\ F_{n+1} & F_n \end{pmatrix}\]
olur. Demek ki \(F_n\), \(Q^n\)’nin birinci satır, ikinci sütundaki elemanıdır; Python’da bu eleman [0, 1] indeksiyle okunur.
import numpy as np
Q = np.array([[1, 1],
[1, 0]])
print(np.linalg.matrix_power(Q, 10))
Q100 = np.linalg.matrix_power(Q, 100)
print(Q100[0, 1]) # int64 taştı!
Q_obj = Q.astype(object) # elemanlar Python tamsayısı olur
print(np.linalg.matrix_power(Q_obj, 100)[0, 1])Çıktı:
[[89 55]
[55 34]]
3736710778780434371
354224848179261915075\(Q^{10}\)’un köşegen dışı elemanından \(F_{10} = 55\) okunur. \(F_{100}\) için ilk sonuç yanlıştır; doğru değer son satırdaki \(354\,224\,848\,179\,261\,915\,075\) sayısıdır. NumPy tamsayı dizilerini varsayılan olarak 64 bitlik int64 türünde tutar. Bu türün en büyük değeri \(2^{63} - 1 \approx 9{,}22 \cdot 10^{18}\)’dir ve \(F_{93}\) bu sınırı aşar. Sınırı aşan bir tamsayı işlemi hata vermez: sonuç \(2^{64}\)’e göre kalanı alınarak sessizce “başa sarar”. astype(object) ile elemanlar, büyüklüğü sınırsız olan Python tamsayılarına çevrilir; hesap yavaşlar ama tam doğru olur. \(\blacksquare\)
Python’ın kendi tamsayıları istediği kadar büyüyebilir (Python ile İlk Adımlar), NumPy’ın int64 dizileri büyüyemez. Ayrıca solve, inv, det ve eig gibi LAPACK’e dayanan fonksiyonlar hesabı her zaman kayan noktalı sayılarla yapar. Tamsayı dizileri önce float64 türüne çevrilir ve sonuç yine yuvarlama hatası taşır: örneğin \(\begin{pmatrix} 2 & 1 \\ 1 & 3 \end{pmatrix}\) tamsayı matrisinin determinantı 5.000000000000001 çıkar. object türündeki diziler ise bu fonksiyonlarda hata verir. Tam aritmetikle matris hesabı gerekiyorsa SymPy ile Sembolik Hesap bölümündeki matrisler kullanılır.
7.2 Lineer Sistemlerin Çözümü
Lineer cebirin merkezindeki problem, \(A\) bir \(n \times n\) matris ve \(b \in \mathbb{R}^n\) iken lineer denklem sistemini \(Ax = b\) biçiminde çözmektir. \(A\) regüler ise çözüm tektir ve NumPy’da tek satırla bulunur: np.linalg.solve(A, b).
solve, elle yaptığımız Gauss yöntemini bir matris ayrışımı olarak düzenler. Satır değişimleriyle birlikte \(A\)’yı bir alt üçgensel \(L\) ve bir üst üçgensel \(U\) matrisinin çarpımına ayırır, sonra iki üçgensel sistemi ileriye ve geriye doğru yerine koyarak çözer. Bu LU ayrışımı yaklaşık \(\tfrac{2}{3}n^3\) aritmetik işlem gerektirir. Her adımda sütundaki mutlak değerce en büyük elemanı köşegene getirmek (kısmi pivotlama, partial pivoting) yuvarlama hatalarının büyümesini önler.
- Katsayılar matrisini
A, sağ tarafıbolarak kayan noktalı diziler biçiminde kurun. x = np.linalg.solve(A, b)ile çözün.- Sağlama yapın:
np.allclose(A @ x, b)ifadesiTrueolmalı, \(r = b - Ax\) artık (residual) vektörü sıfıra yakın olmalıdır.
Örnek 7.2 (Üç Bilinmeyenli Bir Sistem) \[ \begin{aligned} 2x + y - 3z &= 5\\ 3x - 2y + 2z &= 5\\ 5x - 3y - z &= 16 \end{aligned} \]
sistemini solve ile çözünüz ve çözümü sağlayınız.
Çözüm
Üç adımı sırayla uyguluyoruz: katsayıları ve sağ tarafı kuruyor, solve ile çözüyor, sonra sağlama yapıp artığı hesaplıyoruz.
import numpy as np
A = np.array([[2.0, 1.0, -3.0],
[3.0, -2.0, 2.0],
[5.0, -3.0, -1.0]])
b = np.array([5.0, 5.0, 16.0])
x = np.linalg.solve(A, b)
print(x)
print(np.allclose(A @ x, b)) # sağlama
print(b - A @ x) # artık vektörüÇıktı:
[ 1. -3. -2.]
True
[0. 0. 0.]Çözüm \((x, y, z) = (1, -3, -2)\)’dir. Gerçekten \(2 - 3 + 6 = 5\), \(3 + 6 - 4 = 5\) ve \(5 + 9 + 2 = 16\)’dır. Bu sistem Lineer Cebir notlarında Gauss yöntemiyle elle çözülen sistemlerden biridir ve orada da aynı çözüm bulunur. Artığın tam sıfır çıkması bu küçük örneğe özgüdür; genelde \(10^{-16}\) mertebesinde sayılar görürüz. \(\blacksquare\)
solve’un ikinci argümanı bir matris de olabilir. Bu durumda her sütun ayrı bir sağ taraf sayılır ve sonuç, sütunları karşılık gelen çözümler olan bir matristir. LU ayrışımı yalnızca bir kez yapıldığı için bu yol, sistemleri tek tek çözmekten ucuzdur.
import numpy as np
A = np.array([[4.0, 1.0],
[2.0, 3.0]])
B = np.array([[1.0, 0.0, 9.0],
[0.0, 1.0, 7.0]]) # üç sağ taraf, sütun sütun
X = np.linalg.solve(A, B)
print(X)Çıktı:
[[ 0.3 -0.1 2. ]
[-0.2 0.4 1. ]]İlk iki sağ taraf birim matrisin sütunları olduğundan çözümün ilk iki sütunu \(A^{-1}\)’in sütunlarıdır:
\[A^{-1} = \frac{1}{10}\begin{pmatrix} 3 & -1 \\ -2 & 4 \end{pmatrix}.\]
Üçüncü sütun, \((9, 7)\) sağ tarafının çözümü olan \((2, 1)\)’dir.
\(A\) tekil olduğunda solve çoğu zaman LinAlgError hatası verir. Bu hatayı Sınıflar, Hata Yakalama ve Dosyalar bölümündeki try/except yapısıyla yakalayıp sistemi rank ölçütüyle inceleyebiliriz.
Örnek 7.3 (Sonsuz Çözümlü Bir Sistem) \[ \begin{aligned} x + 2y - 3z &= 6\\ 2x - y + 4z &= 2\\ 4x + 3y - 2z &= 14 \end{aligned} \]
sistemini solve ile çözmeyi deneyiniz; başarısız olursa sistemin çözülebilir olup olmadığına rank ölçütüyle karar veriniz.
Çözüm
np.column_stack katsayılar matrisinin yanına \(b\)’yi bir sütun olarak ekler ve genişletilmiş matrisi kurar.
import numpy as np
A = np.array([[1.0, 2.0, -3.0],
[2.0, -1.0, 4.0],
[4.0, 3.0, -2.0]])
b = np.array([6.0, 2.0, 14.0])
try:
x = np.linalg.solve(A, b)
except np.linalg.LinAlgError as err:
print("solve başarısız:", err)
Ab = np.column_stack([A, b]) # genişletilmiş matris [A | b]
print(np.linalg.matrix_rank(A), np.linalg.matrix_rank(Ab))Çıktı:
solve başarısız: Singular matrix
2 2solve hata verdi, çünkü katsayılar matrisi tekildir. matrix_rank hem \(A\)’nın hem de \([A \mid b]\) genişletilmiş matrisinin rankını \(2\) buldu. Çözülebilirliğin rank ölçütüne göre iki rank eşit olduğundan sistem çözülebilirdir. Homojen sistemin çözüm uzayı \(3 - 2 = 1\) boyutludur, yani bir serbest değişken ve sonsuz çok çözüm vardır. Gauss yöntemiyle bulunan genel çözüm \(a \in \mathbb{R}\) olmak üzere \((2 - a,\ 2 + 2a,\ a)\)’dır; örneğin \(a = 1\) için \((1, 4, 1)\) sistemi sağlar: \(1 + 8 - 3 = 6\), \(2 - 4 + 4 = 2\), \(4 + 12 - 2 = 14\). \(\blacksquare\)
7.3 Ters Matris, Determinant ve Rank
Lineer Cebir’de bir kare matrisin tersini, determinantını ve rankını elle hesaplamayı öğrendik. NumPy’daki karşılıkları np.linalg.inv, np.linalg.det ve np.linalg.matrix_rank fonksiyonlarıdır.
import numpy as np
A = np.array([[2.0, 1.0, 1.0],
[1.0, 3.0, 2.0],
[1.0, 0.0, 0.0]])
A_inv = np.linalg.inv(A)
print(A_inv)
print(np.linalg.det(A))
print(np.linalg.matrix_rank(A))
print(np.allclose(A @ A_inv, np.eye(3)))Çıktı:
[[ 0. 0. 1.]
[-2. 1. 3.]
[ 3. -1. -5.]]
-1.0
3
Trueinv ve det de LU ayrışımına dayanır: determinant, \(U\)’nun köşegen elemanlarının çarpımıdır ve satır değişimlerinin sayısı tekse işaret değiştirir. Bu matrisin determinantı \(-1\) olduğu ve elemanları tamsayı olduğu için tersinin elemanları da tamsayıdır (Lineer Cebir, tersin ek matrisle formülü). Elle hesabın aksine bilgisayar her adımda yuvarlama yapar. Bu küçük örnekte görünür bir hata oluşmadı, ama genel durumda A @ A_inv birim matrise yalnızca yaklaşık olarak eşittir; karşılaştırmayı bu yüzden == ile değil np.allclose ile yapıyoruz.
Örnek 7.4 (Sessizce Yanlış Bir Çözüm) \[A = \begin{pmatrix} 0{,}1 & 0{,}2 & 0{,}3 \\ 0{,}4 & 0{,}5 & 0{,}6 \\ 0{,}7 & 0{,}8 & 0{,}9 \end{pmatrix}, \qquad b = \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix}\]
olsun. \(Ax = b\) sisteminin çözümü olmadığını gösterip solve’un bu sisteme verdiği yanıtı inceleyiniz.
Çözüm
\(A\)’nın satırları \(R_1, R_2, R_3\) olsun. \(R_1 - 2R_2 + R_3\) birleşiminin bileşenleri
\[0{,}1 - 0{,}8 + 0{,}7 = 0, \qquad 0{,}2 - 1 + 0{,}8 = 0, \qquad 0{,}3 - 1{,}2 + 0{,}9 = 0\]
olduğundan satırlar lineer bağımlıdır ve \(A\) tekildir. Aynı birleşimi denklemlere uygularsak sol taraf \(0\), sağ taraf \(1 - 2 \cdot 0 + 0 = 1\) olur. \(0 = 1\) çelişkisi sistemin çözümü olmadığını gösterir. Şimdi NumPy’a soralım:
import numpy as np
A = np.array([[0.1, 0.2, 0.3],
[0.4, 0.5, 0.6],
[0.7, 0.8, 0.9]])
b = np.array([1.0, 0.0, 0.0])
x = np.linalg.solve(A, b) # hata vermiyor!
print(x)
print(A @ x) # b'ye hiç benzemiyor
print(np.linalg.det(A))
print(np.linalg.matrix_rank(A))Çıktı:
[-4.50359963e+15 9.00719925e+15 -4.50359963e+15]
[ 0.85 0.2 -0.45]
6.661338147750926e-18
2solve hata vermedi ve bileşenleri \(10^{15}\) mertebesinde olan bir “çözüm” döndürdü. Sağlama yapınca \(Ax\)’in \(b\)’ye hiç benzemediği görülüyor. Determinant da \(0\) değil, \(6{,}7 \cdot 10^{-18}\) çıktı. Sebep, \(0{,}1\), \(0{,}2\), … sayılarının ikili sistemde tam olarak gösterilememesidir (Nümerik Analiz). Bellekteki matris, tekil matrisin çok küçük bir bozulmasıdır ve tam olarak tekil değildir. LU ayrışımı \(0\) olması gereken yerde \(1{,}1 \cdot 10^{-16}\) büyüklüğünde bir pivot bulur ve ona böler. matrix_rank ise doğru yanıtı, \(2\)’yi verdi; bu kararı nasıl verdiğini tekil değer ayrışımı kısmında göreceğiz. \(\blacksquare\)
det(A) == 0 sınaması kayan noktalı sayılarla anlamsızdır. Yukarıdaki tekil matrisin determinantı \(0\) çıkmadı. Tersine, \(0{,}1\, I_{20}\) matrisinin determinantı \(10^{-20}\)’dir, ama bu matrisle sistem çözmek son derece kolaydır. Determinant matrisin ölçeğine bağlıdır ve matrisin tekilliğe ne kadar yakın olduğunu ölçmez. Doğru araçlar matrix_rank ve aşağıda tanımlayacağımız koşul sayısıdır.
7.4 Vektör ve Matris Normları
Bir hesabın ne kadar doğru olduğunu söyleyebilmek için vektörlerin ve matrislerin büyüklüğünü tek bir sayıyla ölçmemiz gerekir. Vektör normlarını Analiz 4 notlarında tanımladık: Öklid normu \(\|x\|_2\), 1-normu \(\|x\|_1 = \sum |x_i|\) ve maksimum normu \(\|x\|_\infty = \max |x_i|\) (\(p\)-normları). Hepsini np.linalg.norm hesaplar; ikinci argüman normun türünü seçer.
import numpy as np
x = np.array([3.0, -4.0, 12.0])
print(np.linalg.norm(x)) # Öklid normu (2-normu)
print(np.linalg.norm(x, 1)) # mutlak değerler toplamı
print(np.linalg.norm(x, np.inf)) # en büyük mutlak değerÇıktı:
13.0
19.0
12.0Gerçekten \(\|x\|_2 = \sqrt{9 + 16 + 144} = 13\), \(\|x\|_1 = 3 + 4 + 12 = 19\) ve \(\|x\|_\infty = 12\)’dir. Matrisler için en kullanışlı normlar, bir vektör normundan türetilenlerdir.
Tanım 7.2 (Doğal Matris Normu) \(\mathbb{R}^n\) ve \(\mathbb{R}^m\) üzerinde aynı türden bir \(\|\cdot\|\) vektör normu verilsin (örneğin ikisinde de Öklid normu). \(m \times n\) bir \(A\) matrisi için (\(x \in \mathbb{R}^n\), \(Ax \in \mathbb{R}^m\))
\[\|A\| = \max_{x \ne 0} \frac{\|Ax\|}{\|x\|} = \max_{\|x\| = 1} \|Ax\|\]
sayısına bu vektör normunun ürettiği doğal matris normu denir. \(\|\cdot\|_1\), \(\|\cdot\|_2\) ve \(\|\cdot\|_\infty\) vektör normlarının ürettiği matris normları sırasıyla \(\|A\|_1\), \(\|A\|_2\) ve \(\|A\|_\infty\) ile gösterilir.
Yani \(\|A\|\), \(A\)’nın bir vektörü en fazla kaç katına uzatabildiğidir. Tanımdan her \(x\) için \(\|Ax\| \le \|A\|\,\|x\|\) eşitsizliği, buradan da \(\|AB\| \le \|A\|\,\|B\|\) eşitsizliği çıkar. Doğal normların yanında, elemanların kareleri toplamından gelen bir norm da sık kullanılır.
Tanım 7.3 (Frobenius Normu) \(m \times n\) bir \(A\) matrisinin Frobenius normu
\[\|A\|_F = \Big(\sum_{i=1}^{m} \sum_{j=1}^{n} a_{ij}^2\Big)^{1/2}\]
sayısıdır.
Yani Frobenius normu, matrisi \(mn\) bileşenli uzun bir vektör gibi düşünüp aldığımız Öklid normudur. Bir vektör normundan türetilmez, ama her \(x\) için \(\|Ax\|_2 \le \|A\|_F\,\|x\|_2\) eşitsizliğini sağlar. 1- ve sonsuz normlarının ise basit formülleri vardır.
Önerme 7.1 (Bir ve Sonsuz Normlarının Formülleri) \(n \times n\) bir \(A\) matrisi için
\[\|A\|_1 = \max_{1 \le j \le n} \sum_{i=1}^{n} |a_{ij}|, \qquad \|A\|_\infty = \max_{1 \le i \le n} \sum_{j=1}^{n} |a_{ij}|\]
dir. Yani \(\|A\|_1\) mutlak değerlerle en büyük sütun toplamı, \(\|A\|_\infty\) en büyük satır toplamıdır.
İspat
Sonsuz normu. \(\|x\|_\infty = 1\) olsun. Her \(j\) için \(|x_j| \le 1\) olduğundan her \(i\) için
\[|(Ax)_i| = \Big|\sum_{j} a_{ij} x_j\Big| \le \sum_{j} |a_{ij}|\,|x_j| \le \sum_{j} |a_{ij}|\]
olur; demek ki \(\|Ax\|_\infty\) en büyük satır toplamını aşmaz. Eşitliği göstermek için en büyük toplamın bulunduğu satır \(k\) olsun ve \(a_{kj} \ge 0\) ise \(x_j = 1\), değilse \(x_j = -1\) seçelim. O zaman \(\|x\|_\infty = 1\) ve \((Ax)_k = \sum_j |a_{kj}|\) olur.
1-normu. Her \(x\) için
\[\|Ax\|_1 = \sum_{i} \Big|\sum_{j} a_{ij} x_j\Big| \le \sum_{j} |x_j| \sum_{i} |a_{ij}| \le \Big(\max_{j} \sum_{i} |a_{ij}|\Big)\, \|x\|_1\]
dir. En büyük sütun toplamı \(k\). sütundaysa \(x = e_k\) için eşitlik sağlanır, çünkü \(Ae_k\), \(A\)’nın \(k\). sütunudur. \(\blacksquare\)
\(\|A\|_2\) için bu kadar basit bir formül yoktur; onu tekil değer ayrışımı kısmında \(A\)’nın en büyük tekil değeri olarak tanıyacağız. np.linalg.norm matrislerde de ikinci argümanla çalışır. Dikkat: ikinci argüman verilmezse matrisin 2-normu değil Frobenius normu hesaplanır.
import numpy as np
A = np.array([[1.0, -2.0],
[3.0, 4.0]])
print(np.linalg.norm(A, 1)) # en büyük sütun toplamı
print(np.linalg.norm(A, np.inf)) # en büyük satır toplamı
print(np.linalg.norm(A, 2)) # 2-normu
print(np.linalg.norm(A)) # varsayılan: Frobenius normuÇıktı:
6.0
7.0
5.116672736016927
5.477225575051661\(A\)’nın mutlak değerlerle sütun toplamları \(1 + 3 = 4\) ve \(2 + 4 = 6\), satır toplamları \(1 + 2 = 3\) ve \(3 + 4 = 7\)’dir. Bu yüzden \(\|A\|_1 = 6\) ve \(\|A\|_\infty = 7\)’dir. Frobenius normu \(\|A\|_F = \sqrt{1 + 4 + 9 + 16}\), yani \(\sqrt{30} \approx 5{,}4772\)’dir.
7.5 Koşul Sayısı
Sessizce yanlış çözüm örneğinde (Örnek 7.4) veride \(10^{-17}\) mertebesindeki yuvarlama hataları, çözümü olmayan bir sistemden bileşenleri \(10^{15}\) mertebesinde olan bir “çözüm” üretti. Bir sistemin verisindeki küçük değişikliklerin çözümü ne kadar değiştirebileceğini ölçen sayı koşul sayısıdır.
Tanım 7.4 (Koşul Sayısı) Regüler bir \(A\) kare matrisinin, bir doğal matris normuna göre koşul sayısı (condition number)
\[\kappa(A) = \|A\|\,\|A^{-1}\|\]
sayısıdır. \(A\) tekil ise \(\kappa(A) = \infty\) kabul edilir. NumPy’da np.linalg.cond(A) 2-normuna göre koşul sayısını, np.linalg.cond(A, 1) ve np.linalg.cond(A, np.inf) öteki iki normdakini verir.
Yani koşul sayısı \(A\) ile \(A^{-1}\)’in büyüklüklerinin çarpımıdır ve hiçbir zaman \(1\)’den küçük olmaz: \(1 = \|I\| = \|AA^{-1}\| \le \|A\|\,\|A^{-1}\|\). Koşul sayısı \(1\)’e yakın matrislere iyi koşullu, çok büyük olanlara kötü koşullu denir. Koşul sayısının neyi ölçtüğünü şu önerme söyler.
Önerme 7.2 (Sağ Taraftaki Hatanın Büyümesi) \(A\) regüler, \(b \ne 0\), \(Ax = b\) ve \(A(x + \delta x) = b + \delta b\) olsun. O zaman
\[\frac{\|\delta x\|}{\|x\|} \le \kappa(A)\, \frac{\|\delta b\|}{\|b\|}\]
dir.
İspat
İki eşitliği taraf tarafa çıkarırsak \(A\,\delta x = \delta b\), yani \(\delta x = A^{-1}\delta b\) olur ve buradan \(\|\delta x\| \le \|A^{-1}\|\,\|\delta b\|\) bulunur. Öte yandan \(\|b\| = \|Ax\| \le \|A\|\,\|x\|\) olduğundan
\[\frac{1}{\|x\|} \le \frac{\|A\|}{\|b\|}\]
dir (\(b \ne 0\) olduğundan \(x \ne 0\)’dır). Bu iki eşitsizliği taraf tarafa çarparsak
\[\frac{\|\delta x\|}{\|x\|} \le \|A\|\,\|A^{-1}\|\,\frac{\|\delta b\|}{\|b\|} = \kappa(A)\,\frac{\|\delta b\|}{\|b\|}\]
elde edilir. \(\blacksquare\)
Yani sağ taraftaki göreli hata çözümde en fazla \(\kappa(A)\) katına çıkabilir; \(A\)’nın kendisindeki küçük hatalar için de benzer bir sınır kanıtlanabilir. Kayan noktalı bir sayının göreli hassasiyeti \(\varepsilon \approx 2{,}2 \cdot 10^{-16}\) olduğundan (Python ile İlk Adımlar) veri zaten bu büyüklükte hatalarla saklanır. Bu yüzden pratik kural şudur: \(\kappa(A) \approx 10^k\) ise çözümün yaklaşık \(16\) anlamlı basamağından \(k\) tanesi kaybedilir.
Örnek 7.5 (Neredeyse Paralel İki Doğru) \(x + y = 2\), \(x - y = 0\) sisteminde ve \(x + y = 2\), \(x + 1{,}1\,y = 2{,}1\) sisteminde ikinci denklemin sağ tarafı \(0{,}1\) artırıldığında çözümlerin ne kadar değiştiğini karşılaştırınız ve sonucu koşul sayısıyla açıklayınız.
Çözüm
İki sistemi ve değiştirilmiş sağ tarafları aynı döngüde çözüyoruz. Her satırda önce özgün, sonra değişmiş sistemin çözümü yazılıyor.
import numpy as np
A_good = np.array([[1.0, 1.0],
[1.0, -1.0]]) # x + y = 2, x - y = 0
A_bad = np.array([[1.0, 1.0],
[1.0, 1.1]]) # x + y = 2, x + 1.1y = 2.1
b_good = np.array([2.0, 0.0])
b_bad = np.array([2.0, 2.1])
db = np.array([0.0, 0.1]) # ikinci bileşene 0.1 ekleniyor
for A, b in [(A_good, b_good), (A_bad, b_bad)]:
x = np.linalg.solve(A, b)
x2 = np.linalg.solve(A, b + db)
print(x, x2)Çıktı:
[1. 1.] [1.05 0.95]
[1. 1.] [0. 2.]İlk sistemde çözüm \((1, 1)\)’den \((1{,}05;\ 0{,}95)\)’e, yani sağ taraftaki değişiklikle aynı mertebede kaydı. İkinci sistemde aynı değişiklik çözümü \((1, 1)\)’den \((0, 2)\)’ye taşıdı. İkinci sistemi koşul sayısıyla inceleyelim:
import numpy as np
A = np.array([[1.0, 1.0],
[1.0, 1.1]])
b = np.array([2.0, 2.1])
b2 = np.array([2.0, 2.2])
x = np.linalg.solve(A, b)
x2 = np.linalg.solve(A, b2)
kappa = np.linalg.cond(A)
rel_b = np.linalg.norm(b2 - b) / np.linalg.norm(b)
rel_x = np.linalg.norm(x2 - x) / np.linalg.norm(x)
print(f"koşul sayısı : {kappa:.4f}")
print(f"b'deki göreli değişim : {rel_b:.4f}")
print(f"x'teki göreli değişim : {rel_x:.4f}")
print(f"önermedeki üst sınır : {kappa * rel_b:.4f}")Çıktı:
koşul sayısı : 42.0762
b'deki göreli değişim : 0.0345
x'teki göreli değişim : 1.0000
önermedeki üst sınır : 1.4509\(b\)’deki \(\%3{,}4\)’lük değişiklik \(x\)’te \(\%100\)’lük bir değişikliğe yol açtı. Büyüme katsayısı \(1 / 0{,}0345 \approx 29\)’dur ve önermedeki üst sınırı (Önerme 7.2), yani \(\kappa(A) \approx 42\) katını aşmaz. İlk sistemin katsayılar matrisinin sütunları dik ve aynı uzunlukta olduğundan koşul sayısı \(1\)’dir. \(\blacksquare\)
Aşağıdaki şekil iki durumu geometrik olarak gösteriyor. İki doğru dik açıyla kesiştiğinde birini biraz kaydırmak kesişim noktasını da biraz kaydırır. Doğrular neredeyse paralel olduğunda ise kesişim noktası çok hassastır: küçük bir kayma onu doğrular boyunca uzağa taşır.
Determinantın bu hassasiyeti ölçmediğini görmek için \(0{,}1\, I_{20}\) matrisine bakalım.
import numpy as np
D = 0.1 * np.eye(20)
print(np.linalg.det(D)) # neredeyse sıfır ...
print(np.linalg.cond(D)) # ... ama koşul sayısı en iyi değerde
print(np.linalg.solve(D, np.ones(20))[:4])Çıktı:
1.0000000000000135e-20
1.0
[10. 10. 10. 10.]Determinant \(10^{-20}\)’dir, ama koşul sayısı mümkün olan en iyi değer olan \(1\)’dir ve sistem tam olarak çözülür. Şimdi koşul sayısı çok hızlı büyüyen ünlü bir matris ailesine bakalım.
Örnek 7.6 (Hilbert Matrisi) Elemanları \(h_{ij} = \dfrac{1}{i + j - 1}\) (\(i, j = 1, \dots, n\)) olan \(n \times n\) matrise Hilbert matrisi denir ve \(H_n\) ile gösterilir. \(n = 5, 6, \dots, 12\) için gerçek çözümü \(x = (1, 1, \dots, 1)\) olan \(H_n x = b\) sistemini kurup solve ile ve inv(H) @ b ile çözünüz; iki yöntemin hatalarını ve artıklarını karşılaştırınız.
Çözüm
Gerçek çözüm birlerden oluşsun diye \(b = H_n \mathbf{1}\) alıyoruz. Matrisi broadcasting ile kuruyoruz (Vektörizasyon ve Broadcasting): i[:, None] + i[None, :] ifadesi \(i + j\) değerlerinden oluşan \(n \times n\) diziyi verir. Hataları göreli olarak ölçüyoruz: hata \(\|\hat x - x\|_2 / \|x\|_2\), artık \(\|H_n \hat x - b\|_2 / \|b\|_2\)’dir; burada \(\hat x\) hesaplanan çözümdür.
import numpy as np
def rel_err(approx, exact):
"""Göreli hata: ||approx - exact|| / ||exact||."""
return np.linalg.norm(approx - exact) / np.linalg.norm(exact)
print(" n koşul hata:solve hata:inv artık:solve artık:inv")
for n in range(5, 13):
i = np.arange(1, n + 1)
H = 1.0 / (i[:, None] + i[None, :] - 1) # Hilbert matrisi
x_true = np.ones(n)
b = H @ x_true
x1 = np.linalg.solve(H, b)
x2 = np.linalg.inv(H) @ b
print(f"{n:2d} {np.linalg.cond(H):8.1e}"
f" {rel_err(x1, x_true):9.1e} {rel_err(x2, x_true):8.1e}"
f" {rel_err(H @ x1, b):10.1e} {rel_err(H @ x2, b):9.1e}")Çıktı:
n koşul hata:solve hata:inv artık:solve artık:inv
5 4.8e+05 1.3e-12 2.7e-11 7.9e-17 2.3e-13
6 1.5e+07 4.1e-11 9.7e-10 1.5e-16 8.1e-12
7 4.8e+08 1.8e-09 1.0e-08 6.5e-17 9.2e-10
8 1.5e+10 6.7e-08 4.8e-07 1.6e-16 2.7e-08
9 4.9e+11 3.0e-06 8.4e-05 7.9e-17 1.9e-06
10 1.6e+13 1.2e-05 2.4e-03 1.6e-16 3.4e-05
11 5.2e+14 1.1e-03 1.0e-01 6.7e-17 4.6e-04
12 1.6e+16 1.3e-01 2.7e+00 1.1e-16 4.6e-02Tablodan üç sonuç okunuyor.
- Koşul sayısı \(n\) bir arttıkça yaklaşık \(33\) katına çıkıyor ve \(n = 12\)’de \(10^{16}\)’ya, yani \(1/\varepsilon\) düzeyine ulaşıyor. Bu noktada koşul sayısının kendisi bile ancak kabaca hesaplanabilir; yüksek duyarlıklı bir hesap gerçek değerin \(1{,}7 \cdot 10^{16}\) olduğunu gösterir.
- İki yöntemin de hatası koşul sayısıyla birlikte büyüyor ve kabaca \(\kappa(H_n)\,\varepsilon\) düzeyinde kalıyor; önermenin (Önerme 7.2) söylediği budur.
invile bulunan çözümün hatası her satırda daha büyük. - Asıl fark artıkta.
solve’un artığı her \(n\) için \(10^{-16}\) civarında kalırkeninv’inki koşul sayısıyla büyüyor; \(n = 12\)’deinvile bulunan çözüm denklemi ancak \(\%5\) civarında bir hatayla sağlıyor. \(\blacksquare\)
7.6 Neden Ters Matrisle Çözmüyoruz?
Matematikte \(Ax = b\)’nin çözümünü \(x = A^{-1}b\) diye yazmak doğaldır ve doğrudur. Kodda bunu np.linalg.inv(A) @ b diye yazmak ise hem daha az doğru hem de daha yavaştır. Hilbert örneğinin (Örnek 7.6) sonuçları aşağıdaki şekilde özetleniyor.
solve ve inv ile bulunan çözümler (dikey eksen logaritmik). Solda göreli hata: ikisi de koşul sayısıyla büyür; kesikli çizgi κ(Hn)·ε kaba tahminidir. Sağda göreli artık ‖Hnx − b‖/‖b‖: solve'un artığı 10−16 düzeyinde kalır, inv'inki koşul sayısıyla birlikte büyür.Artıktaki farkın sebebi şudur. solve’un kullandığı kısmi pivotlamalı LU yöntemi pratikte geriye doğru kararlıdır (backward stable): bulduğu \(\hat x\), \(A\)’ya çok yakın bir \(A + \delta A\) matrisi için (\(\|\delta A\| / \|A\|\) oranı \(\varepsilon\) mertebesinde) sistemin tam çözümüdür. Bu yüzden \(b - A\hat x = \delta A\,\hat x\) artığı her zaman \(\varepsilon\) mertebesindedir. Hata ise bu küçük bozulmanın koşul sayısıyla büyümesinden gelir ve hiçbir yöntem onu ortadan kaldıramaz. inv ile hesaplanan ters matris de yaklaşık olarak doğrudur, ama \(\hat x = \widehat{A^{-1}}\, b\) çarpımı bu güvenceyi taşımaz: artığı koşul sayısıyla büyür.
Hız farkını görmek için rastgele bir \(2000 \times 2000\) sistemi iki yolla çözelim. rng.standard_normal((n, n)) elemanları standart normal dağılımdan, yani ortalaması \(0\) ve standart sapması \(1\) olan normal dağılımdan gelen \(n \times n\) bir dizi üretir. time.perf_counter() saniye cinsinden bir saat değeri verir; her yöntemi beş kez çalıştırıp en kısa süreyi alıyoruz.
import time
import numpy as np
rng = np.random.default_rng(1)
n = 2000
A = rng.standard_normal((n, n))
b = rng.standard_normal(n)
def best_time(f, repeat=5):
"""f'yi birkaç kez çalıştırıp en kısa süreyi saniye cinsinden verir."""
times = []
for _ in range(repeat):
t0 = time.perf_counter()
f()
times.append(time.perf_counter() - t0)
return min(times)
t_solve = best_time(lambda: np.linalg.solve(A, b))
t_inv = best_time(lambda: np.linalg.inv(A) @ b)
print(f"solve : {t_solve:.1g} s")
print(f"inv ve @ : {t_inv:.1g} s")
print(f"oran : yaklaşık {round(t_inv / t_solve)}")Çıktı:
solve : 0.1 s
inv ve @ : 0.4 s
oran : yaklaşık 4Bu süreler bu bilgisayarda, tek iş parçacığıyla ölçüldü; başka bir makinede farklı sayılar çıkar. Oran ise kabaca sabittir: LU ayrışımı yaklaşık \(\tfrac{2}{3}n^3\), tersi hesaplamak yaklaşık \(2n^3\) işlem gerektirir; yani inv yaklaşık üç kat iş yapar. Ölçülen oranın bundan biraz büyük çıkması, süreyi yalnızca işlem sayısının değil, belleğe erişim gibi uygulama ayrıntılarının da belirlemesindendir.
Kodda \(A^{-1}b\) görürseniz onu np.linalg.solve(A, b) olarak yazın. Ters matris yalnızca matrisin kendisi gerçekten gerektiğinde (örneğin elemanları yorumlanacaksa) hesaplanır. Aynı \(A\) ile birçok sistem çözülecekse sağ tarafları bir matrisin sütunlarına koyup solve(A, B) çağırın ya da SciPy’daki scipy.linalg.lu_factor ile ayrışımı bir kez hesaplayıp scipy.linalg.lu_solve ile tekrar tekrar kullanın.
7.7 Özdeğerler ve Özvektörler
Özdeğer ve özvektör kavramlarını Lineer Cebir’den biliyoruz: \(v \ne 0\) ve \(Av = \lambda v\). Elle hesapta önce karakteristik polinomun köklerini buluruz. Sayısal hesapta bu yol kullanılmaz, çünkü bir polinomun kökleri katsayılardaki küçük hatalara çok duyarlıdır. LAPACK bunun yerine, matrisi benzerlik dönüşümleriyle adım adım üçgensel biçime yaklaştıran QR algoritmasını kullanır. NumPy’daki karşılığı np.linalg.eig’dir.
import numpy as np
A = np.array([[2.0, 1.0],
[1.0, 2.0]])
w, V = np.linalg.eig(A)
print(w) # özdeğerler
print(V) # sütunlar: birim özvektörler
print(np.allclose(A @ V, V * w)) # her i için A v_i = w_i v_i
print(np.real_if_close(w)) # sanal kısımlar sıfırsa atÇıktı:
[3.+0.j 1.+0.j]
[[ 0.70710678+0.j -0.70710678+0.j]
[ 0.70710678+0.j 0.70710678+0.j]]
True
[3. 1.]eig iki dizi döndürür: w[i] özdeğerine V[:, i] sütunu karşılık gelir, satırı değil. Özvektörler birim uzunluktadır, işaretleri keyfidir ve özdeğerler belli bir sıraya dizilmez. V * w çarpımı broadcasting ile \(V\)’nin \(i\). sütununu w[i] ile çarpar; böylece bütün \(Av_i = \lambda_i v_i\) eşitlikleri tek satırda sınanır.
Çıktıdaki +0.j kısımları sonucun karmaşık türde (complex128) geldiğini gösteriyor. Eski NumPy sürümleri, sanal kısımların hepsi sıfır olduğunda gerçel türde bir dizi döndürürdü. Kullandığımız NumPy 2.5 sürümünde ise eig sonucu, özdeğerlerin hepsi gerçel olsa bile her zaman karmaşık türde verir; NumPy belgelerindeki eig örnekleri de bunu gösterir. Sanal kısımlar sıfırsa np.real_if_close diziyi gerçel türe çevirir; kodun her sürümde aynı sonucu vermesi için bundan sonra bu adımı ekleyeceğiz.
- Özdeğer
w[i]ise ona ait özvektörV[:, i]sütunudur. np.allclose(A @ V, V * w)ile bütün \(Av_i = \lambda_i v_i\) eşitliklerini sınayın.- Köşegenleştirme için
np.linalg.cond(V)değerine bakın: makul büyüklükteyse \(V\) regülerdir ve \(V^{-1}AV\) köşegendir.
Örnek 7.7 (Bir Matrisin Köşegenleştirilmesi) \[A = \begin{pmatrix} 1 & -3 & 3 \\ 3 & -5 & 3 \\ 6 & -6 & 4 \end{pmatrix}\]
matrisinin özdeğerlerini ve özvektörlerini eig ile bulup \(V^{-1}AV\) matrisinin köşegen olduğunu doğrulayınız.
Çözüm
Üç adımı uyguluyoruz. \(D = V^{-1}AV\) matrisini ters almadan, \(VD = AV\) sistemini solve ile çözerek hesaplıyoruz.
import numpy as np
A = np.array([[1.0, -3.0, 3.0],
[3.0, -5.0, 3.0],
[6.0, -6.0, 4.0]])
w, V = np.linalg.eig(A)
w, V = np.real_if_close(w), np.real_if_close(V)
print(w)
print(np.allclose(A @ V, V * w)) # 2. adım: A v_i = w_i v_i
print(np.linalg.cond(V)) # 3. adım: V regüler mi?
D = np.linalg.solve(V, A @ V) # V^{-1} A V, ters almadan
print(np.allclose(D, np.diag(w)))Çıktı:
[ 4. -2. -2.]
True
4.300625249463308
TrueÖzdeğerler \(4\), \(-2\) ve \(-2\)’dir; \(-2\) iki katlıdır. \(V\)’nin koşul sayısı yaklaşık \(4{,}3\) olduğundan sütunları rahatça bağımsızdır ve \(V^{-1}AV = \operatorname{diag}(4, -2, -2)\) eşitliği sağlanır. Bu matris Lineer Cebir notlarında elle köşegenleştirilmişti. Oradaki \(P\) matrisinin sütunları buradaki \(V\)’nin sütunlarından farklıdır, çünkü iki katlı özdeğerin öz uzayı iki boyutludur ve bu uzayda sonsuz çok taban seçilebilir. \(\blacksquare\)
Örnek 7.8 (Köşegenleştirilemeyen Bir Matris) \[B = \begin{pmatrix} -3 & 1 & -1 \\ -7 & 5 & -1 \\ -6 & 6 & -2 \end{pmatrix}\]
matrisinin karakteristik polinomu bir önceki örnekteki \(A\)’nınkiyle aynıdır, ama \(B\) köşegenleştirilemez. eig’in bu matrise verdiği yanıtı inceleyiniz.
Çözüm
Özdeğerleri, özvektörleri ve \(V\)’nin koşul sayısını yazdıralım:
import numpy as np
B = np.array([[-3.0, 1.0, -1.0],
[-7.0, 5.0, -1.0],
[-6.0, 6.0, -2.0]])
w, V = np.linalg.eig(B)
w, V = np.real_if_close(w), np.real_if_close(V)
print(w)
print(V)
print(np.linalg.cond(V))Çıktı:
[ 4. -1.99999996 -2.00000004]
[[ 8.80956881e-17 -7.07106781e-01 7.07106781e-01]
[-7.07106781e-01 -7.07106781e-01 7.07106781e-01]
[-7.07106781e-01 2.62505234e-08 2.62505208e-08]]
71765775.70757502Tam aritmetikte özdeğerler yine \(4\), \(-2\) ve \(-2\)’dir, ama \(-2\)’nin öz uzayı tek boyutludur (Lineer Cebir). NumPy’ın sonucunda iki şey dikkat çekiyor.
- İki katlı özdeğer \(-1{,}99999996\) ve \(-2{,}00000004\) olarak, yani yalnızca \(8\) basamak doğrulukla bulundu. Köşegenleştirilemeyen bir matrisin katlı özdeğeri, verideki \(\varepsilon\) büyüklüğündeki hataları \(\sqrt{\varepsilon} \approx 10^{-8}\) büyüklüğüne çıkarır.
- \(V\)’nin son iki sütunu neredeyse aynı doğrultudadır: ikisi de \((1, 1, 0)\) doğrultusuna çok yakındır ve biri ötekinin yaklaşık \(-1\) katıdır. Bu yüzden \(V\)’nin koşul sayısı \(7 \cdot 10^7\) civarındadır.
Sayısal hesapta köşegenleştirilemezliği, \(V\)’nin tam olarak tekil çıkmasından değil, koşul sayısının çok büyük olmasından anlarız. \(\blacksquare\)
Gerçel bir matrisin özdeğerleri karmaşık da olabilir. Lineer Cebir notlarındaki \(\begin{pmatrix} 1 & -1 \\ 2 & -1 \end{pmatrix}\) matrisinin karakteristik polinomu \(t^2 + 1\)’dir ve özdeğerleri \(\pm i\)’dir. Python’da sanal birim 1j ile yazılır.
import numpy as np
B = np.array([[1.0, -1.0],
[2.0, -1.0]])
w, V = np.linalg.eig(B)
print(w)
print(V[:, 0])
print(np.allclose(B @ V, V * w))Çıktı:
[-9.71445147e-17+1.j -9.71445147e-17-1.j]
[0.40824829+0.40824829j 0.81649658+0.j ]
TrueGerçel kısımlardaki \(-9{,}7 \cdot 10^{-17}\), tam olarak \(0\) olması gereken bir sayının yuvarlama hatasıdır. \(i\) özdeğerine ait özvektör \(\big(\tfrac{1 + i}{2},\ 1\big)\) vektörünün birim uzunluğa ölçeklenmiş hâlidir: bu vektörün uzunluğunun karesi \(\tfrac{1}{2} + 1 = \tfrac{3}{2}\) olduğundan ölçek çarpanı \(\sqrt{2/3} \approx 0{,}8165\)’tir.
7.8 Simetrik Matrisler
Uygulamada en sık karşılaşılan matrisler simetriktir: kuadratik formların matrisleri, ikinci türev (Hesse) matrisleri, kovaryans matrisleri. Simetrik matrislerin özdeğer problemi çok daha iyi huyludur.
Teorem 7.1 (Simetrik Matrisler İçin Spektral Teorem) \(A\) gerçel ve simetrik (\(A^T = A\)) bir \(n \times n\) matris olsun. O zaman \(A\)’nın bütün özdeğerleri gerçeldir ve sütunları \(A\)’nın özvektörlerinden oluşan, \(Q^TQ = I\) eşitliğini sağlayan bir \(Q\) matrisi vardır:
\[A = Q \Lambda Q^T, \qquad \Lambda = \operatorname{diag}(\lambda_1, \dots, \lambda_n).\]
Yani simetrik bir matris, birbirine dik \(n\) doğrultunun her birini kendi özdeğeri kadar uzatır ya da kısaltır; şekil bunu \(\begin{pmatrix} 2 & 1 \\ 1 & 2 \end{pmatrix}\) matrisi için gösteriyor. Teoremin ispatı iç çarpım uzaylarının konusudur; burada sonuçlarını kullanacağız.
Simetrik matrisler için np.linalg.eigh kullanılır. eigh matrisin yalnızca alt üçgenini okur, eig’den daha hızlıdır, özdeğerleri gerçel türde ve küçükten büyüğe sıralı döndürür ve özvektörlerin dik olmasını garanti eder. Matris simetrik değilse eigh uyarı vermeden yanlış sonuç verir; bu yüzden yalnızca simetrik matrislerde kullanılmalıdır.
Örnek olarak, ikinci türevin sonlu farklarla yaklaşımında ortaya çıkan matrisi alalım. \(x_0 < x_1 < \dots < x_{n+1}\) eşit aralıklı noktalar ve \(h = x_{i+1} - x_i\) ise iç noktalarda
\[f''(x_i) \approx \frac{f(x_{i-1}) - 2f(x_i) + f(x_{i+1})}{h^2}\]
yaklaşımı kullanılır (sonlu fark yaklaşımlarını Sayısal Türev ve İntegral bölümünde inceleyeceğiz). Uçlarda \(f(x_0) = f(x_{n+1}) = 0\) ise bu yaklaşımları \(i = 1, \dots, n\) için alt alta yazmak, iç noktalardaki değerlerin vektörünü \(-L/h^2\) matrisiyle çarpmak demektir. \(L\)’nin köşegeninde \(2\), köşegenin hemen üstünde ve altında \(-1\), öteki yerlerinde \(0\) bulunur. \(n = 3\) için:
import numpy as np
L = np.array([[2.0, -1.0, 0.0],
[-1.0, 2.0, -1.0],
[0.0, -1.0, 2.0]])
w, Q = np.linalg.eigh(L)
print(w) # küçükten büyüğe, gerçel
print(2 - np.sqrt(2), 2 + np.sqrt(2)) # beklenen uç değerler
print(np.allclose(Q.T @ Q, np.eye(3))) # Q ortogonal
print(np.allclose(Q @ np.diag(w) @ Q.T, L)) # L = Q Λ Q^TÇıktı:
[0.58578644 2. 3.41421356]
0.5857864376269049 3.414213562373095
True
TrueBu matrisin özdeğerleri \(2 - \sqrt{2}\), \(2\) ve \(2 + \sqrt{2}\)’dir; eigh bunları sırayla verdi. \(Q^TQ = I\) ve \(L = Q\Lambda Q^T\) eşitlikleri de sağlanıyor.
7.9 En Küçük Kareler
Denklem sayısı bilinmeyen sayısından fazla olan bir \(Ax = b\) sisteminin (\(A\) \(m \times n\) ve \(m > n\)) genellikle çözümü yoktur. Ölçüm verisine bir doğru ya da eğri uydururken karşılaştığımız durum budur. Nümerik Analiz notlarındaki en küçük kareler yöntemi, denklemleri tam sağlamaya çalışmak yerine hataların kareleri toplamını en küçük yapar.
Tanım 7.5 (En Küçük Kareler Çözümü) \(A\) bir \(m \times n\) matris ve \(b \in \mathbb{R}^m\) olsun. \(\|Ax - b\|_2\) değerini en küçük yapan bir \(\hat x \in \mathbb{R}^n\) vektörüne \(Ax = b\) sisteminin en küçük kareler çözümü (least squares solution) denir.
Yani \(\hat x\), \(r = b - A\hat x\) artık vektörünün bileşenlerinin kareleri toplamını en küçük yapar. \(A\)’nın sütunları lineer bağımsızsa (\(\operatorname{rank} A = n\)) \(\hat x\) tektir ve \(A^TA\hat x = A^Tb\) normal denklemlerinin çözümüdür. Doğru uydururken kurulan normal denklemler bunun \(n = 2\) hâlidir. NumPy’da np.linalg.lstsq(A, b) dört şey döndürür: çözüm, artık kareleri toplamı, \(A\)’nın rankı ve \(A\)’nın tekil değerleri.
Örnek 7.9 (Dört Noktaya Doğru Uydurma) \((-1;\ 1)\), \((-0{,}1;\ 1{,}099)\), \((0{,}2;\ 0{,}808)\) ve \((1;\ 1)\) noktalarına en küçük kareler anlamında en iyi uyan \(y = ax + b\) doğrusunu lstsq ile bulunuz.
Çözüm
Her nokta bir \(a x_i + b = y_i\) denklemi verir. Bilinmeyenler \(a\) ve \(b\) olduğundan katsayılar matrisinin sütunları \(x_i\) değerleri ve birlerdir:
\[\begin{pmatrix} -1 & 1 \\ -0{,}1 & 1 \\ 0{,}2 & 1 \\ 1 & 1 \end{pmatrix} \begin{pmatrix} a \\ b \end{pmatrix} \approx \begin{pmatrix} 1 \\ 1{,}099 \\ 0{,}808 \\ 1 \end{pmatrix}.\]
import numpy as np
x = np.array([-1.0, -0.1, 0.2, 1.0])
y = np.array([1.0, 1.099, 0.808, 1.0])
A = np.column_stack([x, np.ones_like(x)]) # sütunlar: x ve 1
coef, rss, rank, sv = np.linalg.lstsq(A, y)
a, b = coef
print(f"a = {a:.9f}, b = {b:.9f}")
print("hata kareleri toplamı:", rss)
print("rank:", rank)
print("tekil değerler:", sv)Çıktı:
a = -0.022454212, b = 0.977311355
hata kareleri toplamı: [0.04347042]
rank: 2
tekil değerler: [2.00127829 1.42999483]Doğru \(y \approx -0{,}022454x + 0{,}977311\)’dir. Bu veri Nümerik Analiz notlarında normal denklemlerle elle çözülmüştü; orada bulunan \(a = -0{,}022454212\) ve \(b = 0{,}977311355\) değerleriyle dokuz basamak uyuşuyor. Hata kareleri toplamı da oradaki \(E \approx 0{,}0435\) değeridir. \(\blacksquare\)
Örnek 7.10 (Gürültülü Veriye Parabol Uydurma) \([0, 4]\) aralığında eşit aralıklı \(15\) noktada \(y = 1 + 2t - 0{,}5t^2\) değerlerine standart sapması \(0{,}3\) olan normal dağılımlı gürültü eklenerek bir veri üretiliyor. Bu veriye \(y = c_0 + c_1 t + c_2 t^2\) parabolünü en küçük kareler yöntemiyle uydurunuz.
Çözüm
rng.normal(0, 0.3, t.size) ortalaması \(0\), standart sapması \(0{,}3\) olan normal dağılımdan t.size tane sayı üretir (Rastgele Sayılar ve Monte Carlo Yöntemleri); tohum sabit olduğu için her çalıştırmada aynı veri çıkar. Model katsayılara göre lineerdir, bu yüzden \(1\), \(t\), \(t^2\) sütunlarından oluşan \(15 \times 3\) bir \(M\) matrisiyle \(Mc \approx y\) sistemini kuruyoruz.
import numpy as np
rng = np.random.default_rng(7)
t = np.linspace(0, 4, 15)
y = 1 + 2 * t - 0.5 * t**2 + rng.normal(0, 0.3, t.size)
M = np.column_stack([np.ones_like(t), t, t**2]) # sütunlar: 1, t, t^2
c, rss, rank, sv = np.linalg.lstsq(M, y)
print("katsayılar:", c)
r = y - M @ c # artıklar
print("artık normu:", np.linalg.norm(r))
print("koşul sayıları:", np.linalg.cond(M), np.linalg.cond(M.T @ M))Çıktı:
katsayılar: [ 0.92977395 2.02839706 -0.50492187]
artık normu: 0.7047186321326601
koşul sayıları: 31.06032228313357 964.7436203321174Uydurulan parabol \(y \approx 0{,}930 + 2{,}028t - 0{,}505t^2\)’dir; gürültüye rağmen katsayılar verinin üretildiği \(1\), \(2\) ve \(-0{,}5\) değerlerine yakın çıktı. Son satır önemli bir gözlem içeriyor: \(M^TM\)’nin koşul sayısı \(M\)’ninkinin karesidir, \(31{,}06^2 \approx 964{,}7\). \(\blacksquare\)
Uydurulan parabol, verinin üretildiği eğri ve artıklar aşağıdaki şekilde görülüyor. Bu tür grafikleri bir sonraki bölümde Matplotlib ile çizeceğiz.
lstsq ile uydurulan parabol y ≈ 0,930 + 2,028t − 0,505t2 ve verinin üretildiği 1 + 2t − 0,5t2 eğrisi (kesikli). Dikey parçalar artıklardır; lstsq bunların kareleri toplamını en küçük yapar.\(A^TA\hat x = A^Tb\) sistemini solve ile çözmek kolay görünür, ama 2-normunda \(\kappa(A^TA) = \kappa(A)^2\) olduğundan kaybedilen basamak sayısı iki katına çıkar (Alıştırma 7.9). lstsq normal denklemleri kurmaz, \(A\)’nın tekil değer ayrışımını kullanır. Polinom uydurmanın hazır araçlarını Polinomlar, İnterpolasyon ve Eğri Uydurma bölümünde göreceğiz.
7.10 Tekil Değer Ayrışımı
Bölüm boyunca birkaç kez tekil değerlere gönderme yaptık: matrix_rank, 2-normu, cond ve lstsq hep onlara dayanır. Şimdi bu ayrışımı tanıyalım.
Tanım 7.6 (Tekil Değer Ayrışımı) \(A\) bir \(m \times n\) gerçel matris ve \(p = \min(m, n)\) olsun. \(U\) (\(m \times m\)) ve \(V\) (\(n \times n\)) matrisleri \(U^TU = I\) ve \(V^TV = I\) eşitliklerini sağlasın; \(\Sigma\) ise köşegeninde \(\sigma_1 \ge \sigma_2 \ge \dots \ge \sigma_p \ge 0\) sayıları, öteki yerlerinde \(0\) bulunan \(m \times n\) bir matris olsun.
\[A = U \Sigma V^T\]
yazılışına \(A\)’nın tekil değer ayrışımı (singular value decomposition, SVD), \(\sigma_i\) sayılarına \(A\)’nın tekil değerleri, \(U\)’nun ve \(V\)’nin \(u_i\) ve \(v_i\) sütunlarına sol ve sağ tekil vektörler denir.
Yani her \(i \le p\) için \(Av_i = \sigma_i u_i\)’dir: \(A\), birbirine dik \(v_i\) doğrultularını birbirine dik \(u_i\) doğrultularına gönderir ve \(i\). doğrultuyu \(\sigma_i\) katına uzatır. Spektral teoremden farkı, giriş ve çıkış için iki ayrı dik taban kullanılmasıdır; bu sayede ayrışım her matris için, kare olmayanlar için bile vardır. NumPy’da np.linalg.svd(A) üç dizi döndürür: U, tekil değerleri büyükten küçüğe içeren tek eksenli s dizisi ve \(V\)’nin kendisi değil, transpozu olan Vt.
import numpy as np
A = np.array([[3.0, 0.0],
[4.0, 5.0]])
U, s, Vt = np.linalg.svd(A)
print(U)
print(s) # tekil değerler, büyükten küçüğe
print(Vt) # dikkat: V değil, V^T
print(np.allclose(U @ np.diag(s) @ Vt, A))
print(np.sqrt(45), np.sqrt(5))Çıktı:
[[-0.31622777 -0.9486833 ]
[-0.9486833 0.31622777]]
[6.70820393 2.23606798]
[[-0.70710678 -0.70710678]
[-0.70710678 0.70710678]]
True
6.708203932499369 2.23606797749979\(A^TA = \begin{pmatrix} 25 & 20 \\ 20 & 25 \end{pmatrix}\) matrisinin özdeğerleri \(45\) ve \(5\)’tir ve tekil değerler bunların karekökleridir: \(\sigma_1 = 3\sqrt5\), \(\sigma_2 = \sqrt5\). Sağ tekil vektörler \(v_1 = -(1, 1)/\sqrt2\) ve \(v_2 = (-1, 1)/\sqrt2\), sol tekil vektörler \(u_1 = -(1, 3)/\sqrt{10}\) ve \(u_2 = (-3, 1)/\sqrt{10}\)’dur. Örneğin \(Av_1 = -(3, 9)/\sqrt2 = \sigma_1 u_1\)’dir. Şekil ayrışımın geometrisini gösteriyor: \(A\) birim çemberi bir elipse götürür ve elipsin yarı eksenleri \(\sigma_1 u_1\) ile \(\sigma_2 u_2\)’dir.
svd'nin döndürdüğü gibidir; işaret seçimi keyfidir.Önerme 7.3 (Tekil Değerlerin Söyledikleri) Her gerçel \(A\) matrisinin bir tekil değer ayrışımı vardır ve tekil değerleri tek türlü belirlidir: \(\sigma_1^2, \dots, \sigma_p^2\) sayıları \(A^TA\)’nın özdeğerleridir (varsa \(A^TA\)’nın öteki özdeğerleri \(0\)’dır). Ayrıca şunlar doğrudur.
(i) \(A\)’nın rankı, sıfırdan farklı tekil değerlerinin sayısıdır.
(ii) \(\|A\|_2 = \sigma_1\)’dir.
(iii) \(A\) kare ve regüler ise \(\kappa(A) = \sigma_1 / \sigma_n\)’dir (2-normunda).
Kare olmayan bir matrisin tersi yoktur, ama (iii)’teki oran onun için de anlamlıdır. Sütunları lineer bağımsız \(m \times n\) bir \(A\) matrisinin (\(m > n\)) 2-normundaki koşul sayısı \(\kappa(A) = \sigma_1 / \sigma_n\) olarak tanımlanır ve np.linalg.cond kare olmayan matrislerde bu oranı hesaplar. En küçük kareler kısmında \(M\) için yazdırdığımız koşul sayısı budur. \(M^TM\)’nin özdeğerleri \(\sigma_i^2\) olduğundan \(\kappa(M^TM) = \sigma_1^2 / \sigma_n^2 = \kappa(M)^2\) eşitliği de buradan çıkar.
Sessizce yanlış çözüm örneğindeki (Örnek 7.4) matrise tekil değerleriyle bakalım. compute_uv=False yalnızca tekil değerleri hesaplatır; np.finfo(float).eps makine epsilonunu verir.
import numpy as np
A = np.array([[0.1, 0.2, 0.3],
[0.4, 0.5, 0.6],
[0.7, 0.8, 0.9]])
s = np.linalg.svd(A, compute_uv=False) # yalnız tekil değerler
print(s)
tol = s.max() * max(A.shape) * np.finfo(float).eps
print(tol, np.sum(s > tol), np.linalg.matrix_rank(A))
print(np.isclose(np.linalg.norm(A, 2), s[0]))
print(s[0] / s[-1], np.linalg.cond(A))Çıktı:
[1.68481034e+00 1.06836951e-01 3.65063904e-17]
1.1223091358001955e-15 2 2
True
4.615110719530146e+16 4.615110719530146e+16Üçüncü tekil değer \(3{,}7 \cdot 10^{-17}\)’dir: tam aritmetikte \(0\) olması gereken bir sayının yuvarlama hatası. matrix_rank, \(\sigma_1 \cdot \max(m, n) \cdot \varepsilon \approx 1{,}1 \cdot 10^{-15}\) eşiğinden büyük tekil değerleri sayar ve rankı \(2\) bulur. 2-normu en büyük tekil değere, koşul sayısı \(\sigma_1/\sigma_3 \approx 4{,}6 \cdot 10^{16}\) oranına eşittir. Bu koşul sayısı \(1/\varepsilon\)’dan büyük olduğundan matris sayısal olarak tekildir ve solve’un ürettiği sayılar anlamsızdır.
Tekil değer ayrışımı \(A\)’yı rankı \(1\) olan matrislerin toplamı olarak da yazar:
\[A = \sigma_1 u_1 v_1^T + \sigma_2 u_2 v_2^T + \dots + \sigma_p u_p v_p^T.\]
Tekil değerler büyükten küçüğe sıralı olduğundan bu toplamın ilk terimleri en önemli kısımdır.
Teorem 7.2 (En İyi Düşük Ranklı Yaklaşım) \(A = U\Sigma V^T\) ve \(k < \operatorname{rank} A\) olsun.
\[A_k = \sum_{i=1}^{k} \sigma_i u_i v_i^T\]
matrisi, rankı en fazla \(k\) olan bütün \(B\) matrisleri arasında \(\|A - B\|_2\) değerini en küçük yapar ve \(\|A - A_k\|_2 = \sigma_{k+1}\)’dir.
Bu sonuç Eckart–Young teoremi olarak bilinir. Yani \(A\)’yı \(k\) terimle yaklaşık yazmanın en iyi yolu ayrışımın ilk \(k\) terimini almaktır ve yapılan hata, atılan ilk tekil değerdir. \(m \times n\) bir matrisin \(A_k\) yaklaşımını saklamak için \(mn\) yerine \(k(m + n + 1)\) sayı yeter.
Örnek 7.11 (Bir Görüntünün Düşük Ranklı Yaklaşımı) \(32 \times 32\) bir matrisin elemanlarını \([0, 1]\) aralığındaki gri tonlar olarak düşünelim (\(0\) zemin rengi, \(1\) en yoğun ton). Bir halka ve onu kesen çapraz bir şeritten oluşan görüntünün rank \(2\), \(6\) ve \(12\) yaklaşımlarını hesaplayıp hataları tekil değerlerle karşılaştırınız.
Çözüm
Izgarayı broadcasting ile kuruyoruz: X sütunun, Y satırın konumunu tutar. Halka ile şeridin kesiştiği hücrelerde toplam \(1\)’i aşar; np.clip(z, 0, 1) bir z dizisinin \(0\)’dan küçük elemanlarını \(0\)’a, \(1\)’den büyük elemanlarını \(1\)’e eşitler. U[:, :k] * s[:k] çarpımı \(U\)’nun ilk \(k\) sütununu karşılık gelen tekil değerlerle çarpar; bunu Vt[:k] ile çarpmak ayrışımın ilk \(k\) teriminin toplamını verir.
import numpy as np
n = 32
g = np.linspace(-1, 1, n)
X, Y = g[None, :], g[:, None] # broadcasting ile ızgara
R = np.sqrt(X**2 + Y**2)
ring = np.exp(-((R - 0.6) / 0.12) ** 2) # bir halka
band = (np.abs(X - Y) < 0.15).astype(float) # çapraz bir şerit
img = np.clip(ring + band, 0, 1) # [0, 1] aralığına kırp
U, s, Vt = np.linalg.svd(img)
print("rank:", np.linalg.matrix_rank(img))
print(np.round(s[:6], 3))
for k in (2, 6, 12):
img_k = (U[:, :k] * s[:k]) @ Vt[:k] # ilk k terimin toplamı
err = np.linalg.norm(img - img_k, 2)
print(k, f"{err:.6f}", f"{s[k]:.6f}", k * (2 * n + 1))Çıktı:
rank: 32
[12.502 5.364 4.514 3.85 3.834 3.251]
2 4.514021 4.514021 130
6 2.872318 2.872318 390
12 1.077171 1.077171 780Keskin kenarlı çapraz şerit tek başına neredeyse tam ranklıdır; bu yüzden görüntünün rankı da tamdır (\(32\)). Her satırda hata, Eckart–Young teoreminin (Teorem 7.2) söylediği gibi atılan ilk tekil değer olan \(\sigma_{k+1}\)’e eşit çıktı. Son sütun yaklaşımı saklamak için gereken sayı adedidir: \(12\) terim için \(780\) sayı, özgün görüntünün \(1024\) sayısından azdır. \(\blacksquare\)
Yaklaşımları görüntü olarak çizdiğimizde rank \(2\)’nin yalnızca kaba bir desen verdiği, rank \(12\)’nin ise halkayı ve şeridi tanınır biçimde yakaladığı görülüyor. Çapraz şerit gibi satırlara ve sütunlara paralel olmayan ayrıntılar çok terim gerektirir.
svd ile bulunan rank 2, 6 ve 12 yaklaşımları. Alttaki sayılar 2-normundaki hatalardır ve sırasıyla σ3, σ7, σ13 tekil değerlerine eşittir. Çizimde 0'ın altına ya da 1'in üstüne taşan değerler 0 ve 1'e kırpıldı.7.11 Alıştırmalar
Aşağıdaki alıştırmaları önce kendiniz çözmeyi deneyin; çözümlerdeki kodlar NumPy 2.5 ile çalıştırılmış ve çıktıları olduğu gibi verilmiştir.
Alıştırma 7.1 (Cramer Kuralının Kodu) Cramer kuralını uygulayan bir cramer(A, b) fonksiyonu yazınız ve Üç Bilinmeyenli Bir Sistem örneğindeki (Örnek 7.2) sistemde sonucunu solve ile karşılaştırınız.
Çözüm
Cramer kuralına göre \(x_j = \det(A_j)/\det(A)\)’dır; \(A_j\), \(A\)’nın \(j\). sütununun yerine \(b\) konarak elde edilir. A.copy() ile kopya almak önemlidir; Aj = A yazsaydık Aj yeni bir dizi değil, aynı dizinin ikinci adı olurdu ve A’nın kendisi değişirdi (NumPy Dizileri).
import numpy as np
def cramer(A, b):
"""Ax = b sistemini Cramer kuralıyla çözer (A kare ve regüler)."""
d = np.linalg.det(A)
x = np.empty(len(b))
for j in range(len(b)):
Aj = A.copy()
Aj[:, j] = b # j. sütunun yerine b konur
x[j] = np.linalg.det(Aj) / d
return x
A = np.array([[2.0, 1.0, -3.0],
[3.0, -2.0, 2.0],
[5.0, -3.0, -1.0]])
b = np.array([5.0, 5.0, 16.0])
print(cramer(A, b))
print(np.linalg.solve(A, b))
print(np.allclose(cramer(A, b), np.linalg.solve(A, b)))Çıktı:
[ 1. -3. -2.]
[ 1. -3. -2.]
Trueİki yöntem aynı sonucu veriyor. Ama Cramer kuralı \(n + 1\) determinant hesaplar; her biri yaklaşık \(\tfrac{2}{3}n^3\) işlem tuttuğundan toplam maliyet \(n^4\) mertebesindedir. Bu yüzden yöntem kuramsal olarak değerlidir, ama büyük sistemlerin çözümünde kullanılmaz. \(\blacksquare\)
Alıştırma 7.2 (Dört Denklemli Bir Sistem) \[ \begin{aligned} x + 2y + 2z &= 2\\ 3x - 2y - z &= 5\\ 2x - 5y + 3z &= -4\\ x + 4y + 6z &= 0 \end{aligned} \]
sisteminin çözülebilir olduğunu rank ölçütüyle gösterip çözümünü lstsq ile bulunuz.
Çözüm
Katsayılar matrisi \(4 \times 3\) olduğundan solve kullanılamaz, çünkü solve kare matris ister. Önce rankları karşılaştırıp sonra lstsq’yu çalıştırıyoruz.
import numpy as np
A = np.array([[1.0, 2.0, 2.0],
[3.0, -2.0, -1.0],
[2.0, -5.0, 3.0],
[1.0, 4.0, 6.0]])
b = np.array([2.0, 5.0, -4.0, 0.0])
Ab = np.column_stack([A, b])
print(np.linalg.matrix_rank(A), np.linalg.matrix_rank(Ab))
x, rss, rank, sv = np.linalg.lstsq(A, b)
print(x)
print(np.linalg.norm(A @ x - b))Çıktı:
3 3
[ 2. 1. -1.]
6.556898606395034e-15\(\operatorname{rank} A = \operatorname{rank}\,[A \mid b] = 3\) olduğundan sistem çözülebilirdir; rank bilinmeyen sayısına eşit olduğundan çözüm tektir. Çözülebilir bir sistemde en küçük kareler çözümü gerçek çözümdür: artık normu \(10^{-15}\) mertebesindedir, yani denklemler yuvarlama hatası dışında tam sağlanır. Çözüm \((x, y, z) = (2, 1, -1)\)’dir: \(2 + 2 - 2 = 2\), \(6 - 2 + 1 = 5\), \(4 - 5 - 3 = -4\) ve \(2 + 4 - 6 = 0\). \(\blacksquare\)
Alıştırma 7.3 (Normlar Arasındaki İki Eşitsizlik) \[A = \begin{pmatrix} 1 & -2 & 3 \\ 0 & 4 & -1 \\ 2 & 1 & 5 \end{pmatrix}\]
matrisinin \(\|A\|_1\), \(\|A\|_\infty\), \(\|A\|_2\) ve \(\|A\|_F\) normlarını hesaplayıp \(\|A\|_2 \le \|A\|_F\) ve \(\|A\|_2 \le \sqrt{\|A\|_1\,\|A\|_\infty}\) eşitsizliklerinin sağlandığını doğrulayınız.
Çözüm
Dört normu hesaplayıp iki eşitsizliği sınıyoruz.
import numpy as np
A = np.array([[1.0, -2.0, 3.0],
[0.0, 4.0, -1.0],
[2.0, 1.0, 5.0]])
n1 = np.linalg.norm(A, 1)
n_inf = np.linalg.norm(A, np.inf)
n2 = np.linalg.norm(A, 2)
n_fro = np.linalg.norm(A, "fro")
print(n1, n_inf, n2, n_fro)
print(n2 <= n_fro, n2 <= np.sqrt(n1 * n_inf))Çıktı:
9.0 8.0 6.399487517369997 7.810249675906654
True TrueMutlak değerlerle sütun toplamları \(3\), \(7\), \(9\) ve satır toplamları \(6\), \(5\), \(8\) olduğundan \(\|A\|_1 = 9\) ve \(\|A\|_\infty = 8\)’dir. \(\|A\|_F = \sqrt{61} \approx 7{,}81\) ve \(\sqrt{9 \cdot 8} = \sqrt{72} \approx 8{,}49\) değerleri \(\|A\|_2 \approx 6{,}40\)’tan büyüktür.
İki eşitsizlik her kare matris için doğrudur. Birincisi tekil değerlerden gelir: Frobenius normu dik matrislerle çarpınca değişmediğinden
\[\|A\|_F^2 = \|\Sigma\|_F^2 = \sigma_1^2 + \dots + \sigma_n^2 \ge \sigma_1^2 = \|A\|_2^2\]
dir. İkincisi için \(\|A\|_2^2\)’nin \(A^TA\)’nın en büyük özdeğeri olduğunu kullanalım. \(Mv = \lambda v\) ve \(v \ne 0\) ise \(|\lambda|\,\|v\|_1 = \|Mv\|_1 \le \|M\|_1\,\|v\|_1\) olduğundan her özdeğerin mutlak değeri \(\|M\|_1\)’i aşmaz. \(M = A^TA\) için bu
\[\|A\|_2^2 \le \|A^TA\|_1 \le \|A^T\|_1\,\|A\|_1 = \|A\|_\infty\,\|A\|_1\]
verir; son adımda \(A^T\)’nin sütunlarının \(A\)’nın satırları olduğunu kullandık. \(\blacksquare\)
Alıştırma 7.4 (Vandermonde Matrislerinin Koşul Sayısı) \([0, 1]\) aralığında eşit aralıklı \(n\) nokta \(x_1, \dots, x_n\) için \((i, j)\) elemanı \(x_i^{\,j-1}\) olan \(n \times n\) Vandermonde matrisinin koşul sayısını \(n = 4, 8, 12, 16, 20\) için hesaplayıp her durumda kaybedilmesi beklenen anlamlı basamak sayısını tahmin ediniz.
Çözüm
np.vander(x, increasing=True) sütunları \(1, x, x^2, \dots\) olan Vandermonde matrisini kurar (NumPy Dizileri). Kaybedilen basamak sayısını \(\log_{10}\kappa\) ile tahmin ediyoruz.
import numpy as np
for n in (4, 8, 12, 16, 20):
x = np.linspace(0, 1, n)
V = np.vander(x, increasing=True) # sütunlar: 1, x, x^2, ...
kappa = np.linalg.cond(V)
print(f"n = {n:2d} koşul = {kappa:8.2e} "
f"kaybedilen basamak ≈ {np.log10(kappa):4.1f}")Çıktı:
n = 4 koşul = 9.89e+01 kaybedilen basamak ≈ 2.0
n = 8 koşul = 2.68e+05 kaybedilen basamak ≈ 5.4
n = 12 koşul = 8.83e+08 kaybedilen basamak ≈ 8.9
n = 16 koşul = 3.12e+12 kaybedilen basamak ≈ 12.5
n = 20 koşul = 1.15e+16 kaybedilen basamak ≈ 16.1Koşul sayısı \(n\) ile üstel olarak büyüyor: \(n = 8\)’de yaklaşık \(5\), \(n = 12\)’de yaklaşık \(9\) basamak kaybedilir ve \(n = 20\)’de koşul sayısı \(10^{16}\)’yı geçer; bu matrisle kurulan bir sistemin çözümünde hiçbir basamağa güvenilemez. Bu matris, verilen noktalardan geçen polinomun katsayılarını bulmak için kurulan sistemin matrisidir. Bu yüzden interpolasyon polinomu pratikte bu sistem çözülerek değil, Lagrange ya da Newton biçimiyle hesaplanır. \(\blacksquare\)
Alıştırma 7.5 (Özdeğerlerle Fibonacci Sayıları) Fibonacci örneğindeki (Örnek 7.1) \(Q\) matrisini \(Q = V\Lambda V^{-1}\) biçiminde köşegenleştirip \(F_n\)’yi \(Q^n = V\Lambda^nV^{-1}\) formülüyle hesaplayınız ve sonucu \(n = 10, 40, 70, 80\) için tam değerle karşılaştırınız.
Çözüm
\(Q\)’nun karakteristik polinomu \(t^2 - t - 1\)’dir; özdeğerleri \(\varphi = \tfrac{1 + \sqrt5}{2}\) ve \(-1/\varphi = \tfrac{1 - \sqrt5}{2}\)’dir. \(\Lambda\) köşegen olduğundan \(\Lambda^n\)’yi hesaplamak köşegen elemanların \(n\). kuvvetini almaktır. \(Q\) simetrik olduğundan eigh kullanıyoruz: döndürdüğü \(V\) ortogonaldir, yani \(V^{-1} = V^T\)’dir ve ters matris hesaplamaya gerek kalmaz. eigh özdeğerleri küçükten büyüğe sıraladığı için önce \(-1/\varphi\) gelir. Tam değerleri dtype=object ile tamsayı matris kuvvetinden alıyoruz.
import numpy as np
Q = np.array([[1.0, 1.0],
[1.0, 0.0]])
w, V = np.linalg.eigh(Q)
print(w) # -1/φ ve φ
Q_obj = np.array([[1, 1], [1, 0]], dtype=object)
for n in (10, 40, 70, 80):
Qn = V @ np.diag(w**n) @ V.T # Q^n = V Λ^n V^T
exact = np.linalg.matrix_power(Q_obj, n)[0, 1]
approx = round(Qn[0, 1])
print(n, approx, exact, approx - exact)Çıktı:
[-0.61803399 1.61803399]
10 55 55 0
40 102334155 102334155 0
70 190392490709135 190392490709135 0
80 23416728348467740 23416728348467685 55\(n = 70\)’e kadar sonuç tam doğru. \(F_{80} = 23\,416\,728\,348\,467\,685\) ise \(17\) basamaklı bir sayıdır; float64 yalnızca \(15\) ile \(16\) arası anlamlı basamak taşıyabildiğinden son basamaklar yanlış çıkar ve fark \(55\) olur. Bu hesabın ifade ettiği formül \(F_n = \big(\varphi^n - (-1/\varphi)^n\big)/\sqrt5\) Binet formülüdür. Matematikte tam olan bu formül, kayan noktalı sayılarla büyük \(n\) için tam sonuç veremez. \(\blacksquare\)
Alıştırma 7.6 (Kuvvet Yöntemi) \[A = \begin{pmatrix} 4 & 1 & 0 \\ 1 & 3 & 1 \\ 0 & 1 & 2 \end{pmatrix}\]
matrisinin mutlak değerce en büyük özdeğerini kuvvet yöntemiyle yaklaşık olarak bulup eigh ile karşılaştırınız: \(x_0 = (1, 1, 1)\) vektöründen başlayıp \(x_{k+1} = Ax_k / \|Ax_k\|_2\) adımlarını tekrarlayınız ve özdeğer tahmini olarak \(\lambda_k = x_k^T A x_k\) Rayleigh bölümünü alınız.
Çözüm
\(A\) simetrik olduğundan özvektörleri dik bir \(q_1, q_2, q_3\) tabanı oluşturur. \(x_0 = c_1q_1 + c_2q_2 + c_3q_3\) ise
\[A^kx_0 = c_1\lambda_1^kq_1 + c_2\lambda_2^kq_2 + c_3\lambda_3^kq_3\]
olur. \(|\lambda_1|\) en büyükse ve \(c_1 \ne 0\) ise ilk terim baskın çıkar, yani \(x_k\) birim vektörü \(\pm q_1\)’e yaklaşır. Her adımda normalleştirmek sayıların taşmasını önler.
import numpy as np
A = np.array([[4.0, 1.0, 0.0],
[1.0, 3.0, 1.0],
[0.0, 1.0, 2.0]])
x = np.array([1.0, 1.0, 1.0])
for k in range(1, 31):
y = A @ x
x = y / np.linalg.norm(y) # her adımda birim vektöre indir
lam = x @ A @ x # Rayleigh bölümü
if k % 5 == 0:
print(f"{k:2d} {lam:.12f}")
w, Q = np.linalg.eigh(A)
print("eigh:", w)Çıktı:
5 4.729619855519
10 4.732025279744
15 4.732050539813
20 4.732050804760
25 4.732050807539
30 4.732050807569
eigh: [1.26794919 3. 4.73205081]En büyük özdeğer \(3 + \sqrt3 \approx 4{,}732050808\)’dir ve \(30\) adımda on iki basamak doğru bulundu. Yakınsama hızını \(|\lambda_2 / \lambda_1| = 3 / 4{,}732 \approx 0{,}634\) oranı belirler: Rayleigh bölümünün hatası her adımda yaklaşık bu oranın karesiyle, \(0{,}40\) ile çarpılır. Tabloda beş adımda hata yaklaşık \(0{,}40^5 \approx 0{,}01\) katına iniyor. İki büyük özdeğer birbirine yakın olsaydı yöntem çok yavaşlardı. \(\blacksquare\)
Alıştırma 7.7 (Bir Markov Zincirinin Daimi Durumu) Geçiş matrisi
\[P = \begin{pmatrix} 0{,}5 & 0{,}5 \\ 0{,}1 & 0{,}9 \end{pmatrix}\]
olan Markov zincirinin daimi durum olasılığını, yani \(\pi P = \pi\) ve \(\pi_1 + \pi_2 = 1\) koşullarını sağlayan \(\pi\) satır vektörünü eig ile bulunuz.
Çözüm
\(\pi P = \pi\) eşitliğinin transpozunu alırsak \(P^T\pi^T = \pi^T\) olur: \(\pi^T\), \(P^T\)’nin \(1\) özdeğerine ait bir özvektörüdür. eig birim uzunlukta bir özvektör döndürür; onu bileşenleri toplamı \(1\) olacak biçimde ölçekliyoruz. np.argmin(np.abs(w - 1)), \(1\)’e en yakın özdeğerin indeksini verir.
import numpy as np
P = np.array([[0.5, 0.5],
[0.1, 0.9]])
w, V = np.linalg.eig(P.T) # πP = π ⟺ P^T π^T = π^T
w, V = np.real_if_close(w), np.real_if_close(V)
print(w)
i = np.argmin(np.abs(w - 1)) # özdeğeri 1 olan sütun
pi = V[:, i] / V[:, i].sum() # bileşenler toplamı 1 olsun
print(pi)
print(np.allclose(pi @ P, pi))Çıktı:
[0.4 1. ]
[0.16666667 0.83333333]
True\(P^T\)’nin özdeğerleri \(0{,}4\) ve \(1\)’dir. Daimi durum \(\pi = (1/6,\ 5/6)\), yani yaklaşık \((0{,}1667;\ 0{,}8333)\)’tür; bu, Raslantı Süreçleri notlarında denklem sistemiyle bulunan sonuçla aynıdır. Zincir uzun vadede zamanın altıda birini birinci, altıda beşini ikinci durumda geçirir. \(\blacksquare\)
Alıştırma 7.8 (Bir Elipsin Asal Eksenleri) \(5x^2 + 4xy + 2y^2 = 1\) elipsinin asal eksen doğrultularını ve yarı eksen uzunluklarını, denklemi \(\mathbf{u} = (x, y)^T\) olmak üzere \(\mathbf{u}^T S\, \mathbf{u} = 1\) biçiminde yazıp \(S\) simetrik matrisinin özdeğer ayrışımıyla bulunuz.
Çözüm
\(\mathbf{u}^TS\,\mathbf{u} = s_{11}x^2 + 2s_{12}xy + s_{22}y^2\) olduğundan karma terimin katsayısı \(4\) köşegen dışına ikiye bölünerek yazılır:
\[S = \begin{pmatrix} 5 & 2 \\ 2 & 2 \end{pmatrix}.\]
Spektral teoremle (Teorem 7.1) \(S = Q\Lambda Q^T\) yazıp \(\mathbf{u}' = Q^T\mathbf{u}\) yeni koordinatlarına geçersek denklem \(\lambda_1 x'^2 + \lambda_2 y'^2 = 1\) olur; yarı eksenler \(1/\sqrt{\lambda_1}\) ve \(1/\sqrt{\lambda_2}\)’dir.
import numpy as np
S = np.array([[5.0, 2.0],
[2.0, 2.0]]) # 5x^2 + 4xy + 2y^2 = [x y] S [x y]^T
w, Q = np.linalg.eigh(S)
print(w)
print(Q)
print(1 / np.sqrt(w)) # yarı eksen uzunlukları
theta = np.degrees(np.arctan(Q[1, 1] / Q[0, 1]))
print(theta) # λ = 6 özvektörünün x ekseniyle açısıÇıktı:
[1. 6.]
[[ 0.4472136 -0.89442719]
[-0.89442719 -0.4472136 ]]
[1. 0.40824829]
26.56505117707799Özdeğerler \(1\) ve \(6\)’dır, dolayısıyla yeni koordinatlarda denklem \(x'^2 + 6y'^2 = 1\)’dir. Büyük yarı eksen \(1\) uzunluğundadır ve \(\lambda = 1\)’in özvektörü \((1, -2)/\sqrt5\) doğrultusundadır. Küçük yarı eksen \(1/\sqrt6 \approx 0{,}408\) uzunluğundadır ve \(\lambda = 6\)’nın özvektörü \((2, 1)/\sqrt5\) doğrultusundadır (eigh işaretleri keyfi seçer). \((2, 1)\) doğrultusunun \(x\) ekseniyle yaptığı açı \(\arctan\tfrac12 \approx 26{,}57^\circ\)’dir. Bu, Analitik Geometri notlarındaki \(\tan 2\theta = 2a_{12}/(a_{11} - a_{22}) = 4/3\) formülünün verdiği dönme açısıdır. \(\blacksquare\)
Alıştırma 7.9 (Normal Denklemlerin Kaybettirdiği Basamaklar) \([0, 1]\) aralığında eşit aralıklı \(30\) noktada \(p(t) = 1 + t + t^2 + \dots + t^9\) polinomunun değerleri veriliyor. Bu veriye dokuzuncu dereceden bir polinomu hem lstsq ile hem de normal denklemleri solve ile çözerek uydurunuz ve bulunan katsayıların en büyük hatalarını karşılaştırınız.
Çözüm
Veri tam olarak bir polinomdan geldiği için doğru katsayıların hepsi \(1\)’dir; hatayı \(\max_k |c_k - 1|\) ile ölçüyoruz. Katsayılar matrisi \(M\), sütunları \(1, t, \dots, t^9\) olan \(30 \times 10\) bir Vandermonde matrisidir.
import numpy as np
t = np.linspace(0, 1, 30)
M = np.vander(t, 10, increasing=True) # sütunlar: 1, t, ..., t^9
c_true = np.ones(10)
y = M @ c_true # tam olarak bir polinomun değerleri
c_lstsq = np.linalg.lstsq(M, y)[0]
c_normal = np.linalg.solve(M.T @ M, M.T @ y)
print(f"koşul(M) = {np.linalg.cond(M):.1e}")
print(f"koşul(M^T M) = {np.linalg.cond(M.T @ M):.1e}")
print(f"lstsq hatası : {np.max(np.abs(c_lstsq - c_true)):.1e}")
print(f"normal denk. hatası : {np.max(np.abs(c_normal - c_true)):.1e}")Çıktı:
koşul(M) = 3.5e+06
koşul(M^T M) = 1.2e+13
lstsq hatası : 9.2e-11
normal denk. hatası : 9.3e-04\(M\)’nin koşul sayısı \(3{,}5 \cdot 10^6\), \(M^TM\)’ninki onun karesi olan \(1{,}2 \cdot 10^{13}\)’tür. lstsq katsayıları \(10^{-10}\) mertebesinde bir hatayla bulurken normal denklemlerin hatası \(10^{-3}\) mertebesine çıktı: aşağı yukarı yedi basamak fazladan kaybedildi. \(\blacksquare\)
Alıştırma 7.10 (Tekil Değer Ayrışımıyla Çekirdek) \[M = \begin{pmatrix} 1 & 2 & 3 \\ 4 & 5 & 6 \\ 7 & 8 & 9 \end{pmatrix}\]
matrisinin çekirdeğini, yani \(Mv = 0\) denkleminin çözüm uzayını tekil değer ayrışımı yardımıyla bulunuz.
Çözüm
\(M = U\Sigma V^T\) ise \(Mv_i = \sigma_i u_i\)’dir. Bu yüzden sıfır olan (ya da sayısal olarak sıfır sayılan) tekil değerlere ait sağ tekil vektörler çekirdeği gerer. Tekil değerler büyükten küçüğe sıralı olduğundan aranan vektör Vt’nin son satırıdır.
import numpy as np
M = np.array([[1.0, 2.0, 3.0],
[4.0, 5.0, 6.0],
[7.0, 8.0, 9.0]])
U, s, Vt = np.linalg.svd(M)
print(s)
v = Vt[-1] # en küçük tekil değere ait sağ tekil vektör
print(v / v[0]) # ilk bileşeni 1 olacak biçimde ölçekle
print(np.linalg.norm(M @ v))Çıktı:
[1.68481034e+01 1.06836951e+00 3.33475287e-16]
[ 1. -2. 1.]
5.20740757162067e-16Üçüncü tekil değer \(3{,}3 \cdot 10^{-16}\), yani sayısal olarak sıfırdır; öteki ikisi açıkça sıfırdan farklı olduğundan \(\operatorname{rank} M = 2\) ve çekirdek bir boyutludur. Son sağ tekil vektör, ilk bileşeni \(1\) olacak biçimde ölçeklenince \((1, -2, 1)\) çıkıyor. Gerçekten \(1 - 4 + 3 = 0\), \(4 - 10 + 6 = 0\) ve \(7 - 16 + 9 = 0\)’dır; çekirdek \(\{t\,(1, -2, 1) : t \in \mathbb{R}\}\) doğrusudur. \(\blacksquare\)
Alıştırma 7.11 (En İyi Rank Bir Yaklaşım) \[A = \begin{pmatrix} 3 & 2 & 2 \\ 2 & 3 & -2 \end{pmatrix}\]
matrisine 2-normunda en yakın olan, rankı \(1\) olan \(A_1\) matrisini bulup \(\|A - A_1\|_2\) değerini hesaplayınız.
Çözüm
Eckart–Young teoremine (Teorem 7.2) göre aranan matris \(A_1 = \sigma_1 u_1 v_1^T\)’dir ve hata \(\sigma_2\)’dir. np.outer(u, v) iki vektörün dış çarpımı olan \(uv^T\) matrisini verir.
import numpy as np
A = np.array([[3.0, 2.0, 2.0],
[2.0, 3.0, -2.0]])
U, s, Vt = np.linalg.svd(A)
print(s)
A1 = s[0] * np.outer(U[:, 0], Vt[0]) # σ1 u1 v1^T
print(np.round(A1, 12))
print(np.linalg.norm(A - A1, 2))Çıktı:
[5. 3.]
[[ 2.5 2.5 -0. ]
[ 2.5 2.5 -0. ]]
2.9999999999999996Tekil değerler \(5\) ve \(3\)’tür; gerçekten \(AA^T = \begin{pmatrix} 17 & 8 \\ 8 & 17 \end{pmatrix}\) matrisinin özdeğerleri \(25\) ve \(9\)’dur. \(u_1 = (1, 1)/\sqrt2\) ve \(v_1 = (1, 1, 0)/\sqrt2\) olduğundan (işaretler birlikte değişebilir)
\[A_1 = 5 \cdot \frac{1}{2}\begin{pmatrix} 1 & 1 & 0 \\ 1 & 1 & 0 \end{pmatrix} = \begin{pmatrix} 2{,}5 & 2{,}5 & 0 \\ 2{,}5 & 2{,}5 & 0 \end{pmatrix}\]
bulunur. Çıktıdaki -0., \(-10^{-16}\) gibi çok küçük negatif bir sayının yuvarlanmasından doğan işaretli sıfırdır; IEEE 754’te \(-0\) ile \(0\) eşittir. Hata \(\|A - A_1\|_2 = 3 = \sigma_2\)’dir; çıktıdaki \(2{,}9999999999999996\) yine yuvarlama hatasıdır. \(\blacksquare\)
Bu bölümde NumPy’ın lineer cebir araçlarını ve onları güvenle kullanmanın kurallarını gördük: sistemleri inv ile değil solve ile çözmek, tekilliği determinantla değil rank ve koşul sayısıyla sınamak, simetrik matrislerde eigh’i seçmek, aşırı belirli sistemlerde normal denklemler yerine lstsq’ya başvurmak. Sonuçları bu bölümde sık sık şekillerle gösterdik; bir sonraki bölümde, Matplotlib ile Grafik Çizimi bölümünde, bu tür grafikleri kendimiz çizmeyi öğreneceğiz.