6 Vektörizasyon ve Broadcasting
Matematikte bir fonksiyonu bir sayı kümesinin her elemanına uygularken “her \(i\) için \(y_i = x_i^2 - \sin x_i\)” deyip geçeriz; kimse bunu tek tek yazmaz. Kontrol Yapıları ve Fonksiyonlar bölümünde aynı işi bir for döngüsüyle yapmıştık. NumPy Dizileri bölümünde ise dizileri kurmayı, şekillerini değiştirmeyi ve indekslemeyi öğrendik. Bu bölümde dizilerle döngü yazmadan hesap yapacağız:
import math
import numpy as np
# Döngüyle: her eleman ayrı ayrı işlenir
xs = [0.0, 0.5, 1.0, 1.5, 2.0]
ys = []
for v in xs:
ys.append(v**2 - math.sin(v))
print([round(v, 8) for v in ys]) # NumPy gibi 8 basamakla yazdır
# Vektörel: bütün dizi tek ifadeyle işlenir
x = np.array(xs)
y = x**2 - np.sin(x)
print(y)Çıktı:
[0.0, -0.22942554, 0.15852902, 1.25250501, 3.09070257]
[ 0. -0.22942554 0.15852902 1.25250501 3.09070257]İki hesap aynı sayıları verir, ama ikincisi matematikteki yazılışa çok daha yakındır. Bir dizinin bütün elemanlarını tek bir ifadeyle işleme tekniğine vektörizasyon (vectorization) denir. Vektörel kod hem kısadır hem de çok hızlıdır, çünkü elemanlar üzerindeki döngü Python yorumlayıcısında değil, NumPy’ın derlenmiş C kodunda döner. Bölümün ikinci büyük konusu broadcasting kurallarıdır: şekilleri farklı iki diziyi, örneğin bir matrisle bir satırı, kopya oluşturmadan birlikte işlemeyi sağlarlar.
6.1 Evrensel Fonksiyonlar
Vektörel hesabın yapı taşı, bir diziye eleman eleman uygulanan fonksiyonlardır.
import numpy as np
x = np.linspace(0, np.pi, 5)
print(x)
print(np.sin(x))
a = np.array([1, 2, 3, 4])
b = np.array([10, 20, 30, 40])
print(a + b, np.add(a, b)) # aynı işlem, iki yazılış
print(a * b, a**2 / b)
print(np.maximum(a, 3)) # her i için max(a[i], 3)
print(np.sqrt(a))Çıktı:
[0. 0.78539816 1.57079633 2.35619449 3.14159265]
[0.00000000e+00 7.07106781e-01 1.00000000e+00 7.07106781e-01
1.22464680e-16]
[11 22 33 44] [11 22 33 44]
[ 10 40 90 160] [0.1 0.2 0.3 0.4]
[3 3 3 4]
[1. 1.41421356 1.73205081 2. ]np.sin(x) beş elemanın beşinin de sinüsünü verdi. Son değer \(\sin \pi = 0\) olması gerekirken \(1{,}22 \cdot 10^{-16}\) çıktı: np.pi sayısı \(\pi\)’nin kendisi değil, ona en yakın makine sayısıdır (bkz. Python ile İlk Adımlar). Bu küçük sayıyı da gösterebilmek için NumPy ikinci satırdaki bütün diziyi bilimsel gösterimle yazdı. a + b ile np.add(a, b) aynı işlemdir; +, *, **, / işleçleri NumPy dizilerinde arka planda bu fonksiyonları çağırır.
Tanım 6.1 (Evrensel Fonksiyon) Aynı şekle sahip dizilere eleman eleman uygulanan ve sonucu yine aynı şekle sahip bir dizi olarak veren NumPy fonksiyonuna evrensel fonksiyon (universal function, kısaca ufunc) denir. Tek girdili bir ufunc \(f\) için \(y = f(x)\) dizisinin elemanları \(y_i = f(x_i)\), iki girdili bir \(g\) için \(z = g(x, y)\) dizisinin elemanları \(z_i = g(x_i, y_i)\) olur.
Yani bir ufunc, sayılar için tanımlı bir fonksiyonu dizilere “kendiliğinden” genişletir; döngüyü biz değil, NumPy kurar. Girdi tek bir sayı da olabilir; o zaman sonuç da tek bir sayıdır. En sık kullanılanlar şunlardır:
| Tür | Evrensel fonksiyonlar |
|---|---|
| Aritmetik | np.add (+), np.subtract (-), np.multiply (*), np.divide (/), np.floor_divide (//), np.mod (%), np.power (**) |
| Üstel ve logaritma | np.exp, np.log, np.log2, np.log10, np.sqrt |
| Trigonometrik | np.sin, np.cos, np.tan, np.arcsin, np.arctan, np.arctan2, np.hypot, np.sinh, np.cosh |
| Yuvarlama ve işaret | np.abs, np.sign, np.floor, np.ceil, np.rint (en yakın tamsayıya) |
| Karşılaştırma | np.maximum, np.minimum, np.greater (>), np.less (<), np.equal (==) |
Sık kullanılan np.round (verilen basamağa yuvarlama) ve np.isclose (yaklaşık eşitlik) da eleman eleman çalışır, ama ufunc değildir. Bu yüzden ufunc’ların birazdan göreceğimiz reduce, accumulate gibi yöntemleri onlarda yoktur.
Python ile İlk Adımlar bölümündeki math modülünün fonksiyonları ise yalnız tek bir sayı kabul eder:
import math
import numpy as np
x = np.linspace(0, 1, 3)
try:
math.sin(x) # math yalnız tek sayı kabul eder
except TypeError as err:
print("Hata:", err)
print(np.sin(x)) # np.sin dizinin her elemanına uygulanır
print(math.sin(0.5), np.sin(0.5))Çıktı:
Hata: only 0-dimensional arrays can be converted to Python scalars
[0. 0.47942554 0.84147098]
0.479425538604203 0.479425538604203Kural basittir: dizilerle çalışırken matematik fonksiyonlarını np içinden alın. Tek bir sayı için math.sin ile np.sin aynı değeri verir.
Python’da math.sqrt(-1) ya da 1 / 0 hata verip programı durdurur. NumPy ise dizinin geri kalanını hesaplamaya devam eder. Geçersiz sonucun yerine nan (“not a number”, sayı değil), sıfıra bölmede ya da taşmada inf veya -inf yazar ve yalnız bir RuntimeWarning uyarısı basar. Aşağıda np.errstate bu uyarıları susturuyor ki çıktı yalnız sonuçları göstersin:
import numpy as np
x = np.array([4.0, 1.0, 0.0, -1.0])
with np.errstate(divide="ignore", invalid="ignore"):
r = np.sqrt(x)
lg = np.log(x)
inv = 1 / x
print(r)
print(lg)
print(inv)
print(np.isnan(r), np.isfinite(lg))Çıktı:
[ 2. 1. 0. nan]
[1.38629436 0. -inf nan]
[ 0.25 1. inf -1. ]
[False False False True] [ True True False False]nan sonraki her işleme bulaşır: nan + 1 de, nan içeren bir dizinin toplamı da nan’dır. Şüpheli sonuçları np.isnan ve np.isfinite ile ayıklayın; x == np.nan her zaman False verir.
İki girdili ufunc’ların birkaç yararlı yöntemi de vardır. reduce işlemi bütün diziye art arda uygulayıp tek sayıya indirir, accumulate ara sonuçları da tutar, outer ise bütün \((x_i, y_j)\) çiftleri için tablo kurar:
import numpy as np
n = np.arange(1, 8)
print(np.add.reduce(n)) # 1 + 2 + ... + 7
print(np.add.accumulate(n)) # kısmi toplamlar
print(np.multiply.accumulate(n)) # 1!, 2!, ..., 7!
print(np.multiply.outer(n[:4], n[:4]))Çıktı:
28
[ 1 3 6 10 15 21 28]
[ 1 2 6 24 120 720 5040]
[[ 1 2 3 4]
[ 2 4 6 8]
[ 3 6 9 12]
[ 4 8 12 16]]np.multiply.accumulate ilk yedi faktöriyeli, np.multiply.outer ise küçük bir çarpım tablosunu tek satırda verdi. Uygulamada np.add.reduce yerine eşdeğeri olan np.sum, np.add.accumulate yerine np.cumsum yazılır; bunları birazdan indirgeme işlemleriyle birlikte göreceğiz.
Örnek 6.1 (Kartezyenden Kutupsal Koordinatlara) \((1, 1)\), \((0, 2)\), \((-1, 0)\), \((-1, -1)\) ve \((3, -4)\) noktalarının kutupsal koordinatlarını, yani \(r = \sqrt{x^2 + y^2}\) uzaklığını ve \((-\pi, \pi]\) aralığındaki \(\theta\) açısını döngü kullanmadan hesaplayın.
Çözüm
Apsisleri ve ordinatları iki ayrı dizide tutarız. \(r\) için np.hypot(x, y) evrensel fonksiyonu \(\sqrt{x_i^2 + y_i^2}\) değerini her nokta için hesaplar. Açı için np.arctan(y / x) kullanamayız. \(x = 0\) olan \((0, 2)\) noktasında sıfıra bölmüş oluruz. Ayrıca \(\arctan\) yalnız \((-\pi/2, \pi/2)\) aralığında değer aldığından \((-1, -1)\) ile \((1, 1)\) noktalarına aynı açıyı verir. np.arctan2(y, x) ise noktanın hangi bölgede olduğuna da bakarak doğru açıyı verir.
import numpy as np
x = np.array([1.0, 0.0, -1.0, -1.0, 3.0])
y = np.array([1.0, 2.0, 0.0, -1.0, -4.0])
r = np.hypot(x, y) # sqrt(x**2 + y**2)
theta = np.arctan2(y, x) # (-pi, pi] aralığında açı
print(r)
print(np.round(np.degrees(theta), 2))
# Geri dönüşüm: x = r cos(theta), y = r sin(theta)
print(np.allclose(r * np.cos(theta), x),
np.allclose(r * np.sin(theta), y))Çıktı:
[1.41421356 2. 1. 1.41421356 5. ]
[ 45. 90. 180. -135. -53.13]
True TrueAçılar derece cinsinden \(45°\), \(90°\), \(180°\), \(-135°\) ve yaklaşık \(-53{,}13°\)’dir. \((-1, -1)\) noktası üçüncü bölgede olduğu için açısı \(-135°\) çıktı, \((1, 1)\) ile karışmadı. Son satır \(x = r\cos\theta\) ve \(y = r\sin\theta\) ile geri dönüşümün bütün noktalar için tuttuğunu np.allclose ile doğruluyor. Bu fonksiyon iki dizinin eleman eleman yuvarlama hatası kadar yakın olup olmadığını sınar. \(\blacksquare\)
6.2 İndirgeme İşlemleri ve axis Parametresi
Bir dizinin toplamı, en büyük elemanı ya da ortalaması gibi sayılar, diziyi tek bir sayıya indirgeyen işlemlerle bulunur.
import numpy as np
v = np.array([3.0, -1.0, 4.0, 1.0, -5.0, 9.0])
print(v.sum(), v.prod(), v.mean())
print(v.min(), v.max(), v.argmin(), v.argmax())
print(v.std(), np.median(v))
print(np.cumsum(v))Çıktı:
11.0 540.0 1.8333333333333333
-5.0 9.0 4 5
4.336537277085896 2.0
[ 3. 2. 6. 7. 2. 11.]argmin ve argmax en küçük ve en büyük elemanın değerini değil, indeksini verir: en büyük eleman olan \(9\), dizinin 5 numaralı yerindedir. cumsum kısmi toplamları, std standart sapmayı verir. std varsayılan olarak \(n\)’ye böler; örneklem standart sapması için v.std(ddof=1) yazılır. Tamsayı dizilerinde prod kolayca int64 sınırını aşar. NumPy Dizileri bölümünde gördüğümüz gibi bu taşma sessizdir ve örneğin \(21!\) negatif bir sayı olarak çıkar.
Çok boyutlu dizilerde çoğu zaman bütün elemanların toplamını değil, satır ya da sütun toplamlarını isteriz. Hangi eksen boyunca toplanacağını, yani NumPy Dizileri bölümünde tanıdığımız eksenlerden hangisinin kullanılacağını axis parametresi belirler.
Tanım 6.2 (Bir Eksen Boyunca İndirgeme) Şekli \((n_0, n_1, \ldots, n_{d-1})\) olan bir \(A\) dizisinin \(k\) numaralı eksen boyunca toplamı, A.sum(axis=k), \(k\)-ıncı indeks üzerinden toplanarak elde edilen dizidir. Bu dizinin şekli, \(A\)’nın şeklinden \(n_k\) çıkarılarak bulunur. Bir matris (\(d = 2\)) için A.sum(axis=0) dizisinin \(j\)’inci elemanı \(\sum_{i} A_{ij}\), A.sum(axis=1) dizisinin \(i\)’inci elemanı ise \(\sum_{j} A_{ij}\) olur. prod, mean, std, min, max, argmin ve argmax da axis parametresini aynı anlamda kullanır. axis verilmezse bütün elemanlar üzerinden tek bir sonuç hesaplanır.
Yani axis ile adı verilen eksen sonuçtan kaybolur. Bir matriste axis=0 satır indeksi \(i\)’yi yok eder: her sütun aşağı doğru toplanır ve sütun sayısı kadar sonuç çıkar. axis=1 ise sütun indeksi \(j\)’yi yok eder ve her satırın toplamını verir.
axis=0: toplama 0. eksen boyunca, yani satırlar üzerinden yapılır, her sütun tek sayıya iner (ilk sütun 0 + 4 + 8 = 12) ve sonuç (4,) şeklindedir. Sağda axis=1: her satır tek sayıya iner (ilk satır 0 + 1 + 2 + 3 = 6) ve sonuç (3,) şeklindedir. Adı verilen eksen sonuçtan kaybolur.import numpy as np
A = np.arange(12).reshape(3, 4)
print(A)
print(A.sum(axis=0)) # sütun toplamları, şekil (4,)
print(A.sum(axis=1)) # satır toplamları, şekil (3,)
print(A.sum()) # bütün elemanların toplamı
print(A.max(axis=0), A.argmax(axis=1))
print(A.sum(axis=1, keepdims=True))
print(A.cumsum(axis=1)) # kısmi toplamlar, şekil (3, 4) korunur
B = np.arange(24).reshape(2, 3, 4)
print(B.sum(axis=1).shape, B.sum(axis=(0, 2)))Çıktı:
[[ 0 1 2 3]
[ 4 5 6 7]
[ 8 9 10 11]]
[12 15 18 21]
[ 6 22 38]
66
[ 8 9 10 11] [3 3 3]
[[ 6]
[22]
[38]]
[[ 0 1 3 6]
[ 4 9 15 22]
[ 8 17 27 38]]
(2, 4) [ 60 92 124]keepdims=True verilirse kaybolan eksen silinmez, uzunluğu 1 olarak kalır: satır toplamları (3,) yerine (3, 1) şeklinde bir sütun olarak gelir. Bu, birazdan göreceğimiz broadcasting ile birlikte çok işe yarar. cumsum ve cumprod ise indirgeme değildir: axis ile verilen eksen boyunca kısmi toplamları (çarpımları) hesaplarlar ve dizinin şeklini korurlar. Yukarıda A.cumsum(axis=1) her satırın kısmi toplamlarını verdi ve sonuç yine (3, 4) şeklinde çıktı. axis verilmezse dizi önce düzleştirilir; A.cumsum() 12 elemanlı, tek eksenli bir dizidir. Son satırda üç boyutlu bir dizi var. axis=1 ortadaki uzunluğu 3 olan ekseni yok ettiği için sonuç (2, 4) şeklindedir. axis=(0, 2) ise iki ekseni birden yok eder.
Örnek 6.2 (Dürer’in Sihirli Karesi) Albrecht Dürer’in 1514 tarihli Melencolia I gravüründeki
\[ M = \begin{pmatrix} 16 & 3 & 2 & 13 \\ 5 & 10 & 11 & 8 \\ 9 & 6 & 7 & 12 \\ 4 & 15 & 14 & 1 \end{pmatrix} \]
karesinde bütün satır ve sütun toplamlarının, iki köşegenin toplamının ve dört \(2 \times 2\) köşe bloğunun toplamlarının \(34\) olduğunu NumPy ile gösterin.
Çözüm
Satır toplamları axis=1, sütun toplamları axis=0 ile gelir. Köşegen toplamı izdir (np.trace). Ters köşegen için matrisi np.fliplr ile soldan sağa çevirip izini alırız. Dört köşe bloğu için \(4 \times 4\) matrisi reshape(2, 2, 2, 2) ile dört indeksli bir diziye çeviririz. Yeni dizide \(M_{ij}\) elemanı, \(i = 2a + b\) ve \(j = 2c + e\) olmak üzere \((a, b, c, e)\) konumundadır. \((a, c)\) hangi blokta olduğumuzu, \((b, e)\) blok içindeki konumu gösterir. Blok içi indeksler 1 ve 3 numaralı eksenlerdir, bu yüzden onlar üzerinden toplarız.
import numpy as np
M = np.array([[16, 3, 2, 13],
[5, 10, 11, 8],
[9, 6, 7, 12],
[4, 15, 14, 1]])
print("satırlar :", M.sum(axis=1))
print("sütunlar :", M.sum(axis=0))
print("köşegen :", np.trace(M))
print("ters köş. :", np.trace(np.fliplr(M)))
print("çeyrekler :", M.reshape(2, 2, 2, 2).sum(axis=(1, 3)).ravel())Çıktı:
satırlar : [34 34 34 34]
sütunlar : [34 34 34 34]
köşegen : 34
ters köş. : 34
çeyrekler : [34 34 34 34]Hepsi \(34\)’tür. Bu sayı rastgele değildir: \(1, 2, \ldots, 16\) sayılarının toplamı \(136\)’dır. Dört satırın toplamı eşitse ortak değer zorunlu olarak \(136 / 4 = 34\) olur. Yani \(1\)’den \(16\)’ya kadar sayılarla kurulan her \(4 \times 4\) sihirli karenin sabiti \(34\)’tür. \(\blacksquare\)
Örnek 6.3 (Basel Serisinin Kısmi Toplamları) \(S_n = \sum_{k=1}^{n} 1/k^2\) kısmi toplamlarının hepsini \(n = 10^6\)’ya kadar tek bir cumsum çağrısıyla hesaplayın. \(n = 10, 100, \ldots, 10^6\) için \(\pi^2/6 - S_n\) farkını yazdırıp farkın nasıl küçüldüğünü yorumlayın.
Çözüm
Kontrol Yapıları ve Fonksiyonlar bölümünde bu toplamları iç içe döngülerle, her \(n\) için baştan hesaplamıştık. Şimdi hepsini tek seferde buluruz. \(k = 1, 2, \ldots, 10^6\) dizisini np.arange ile kurarız. 1.0 / k**2 dizisi serinin terimlerini, bunun cumsum’ı da bütün kısmi toplamları verir. Python indeksleri 0’dan başladığı için \(S_n\) değeri S[n - 1] konumundadır.
import numpy as np
N = 10**6
k = np.arange(1, N + 1)
S = np.cumsum(1.0 / k**2) # S[n - 1] = 1 + 1/4 + ... + 1/n^2
target = np.pi**2 / 6
for n in [10, 100, 1000, 10**4, 10**5, 10**6]:
gap = target - S[n - 1]
print(f"n = {n:>7} S_n = {S[n - 1]:.10f} "
f"fark = {gap:.3e} n * fark = {n * gap:.4f}")Çıktı:
n = 10 S_n = 1.5497677312 fark = 9.517e-02 n * fark = 0.9517
n = 100 S_n = 1.6349839002 fark = 9.950e-03 n * fark = 0.9950
n = 1000 S_n = 1.6439345667 fark = 9.995e-04 n * fark = 0.9995
n = 10000 S_n = 1.6448340718 fark = 1.000e-04 n * fark = 1.0000
n = 100000 S_n = 1.6449240669 fark = 1.000e-05 n * fark = 1.0000
n = 1000000 S_n = 1.6449330668 fark = 1.000e-06 n * fark = 1.0000Son sütundaki \(n \cdot (\pi^2/6 - S_n)\) değerleri \(1\)’e gidiyor, yani kalan yaklaşık \(1/n\)’dir. Bu, integral testinden gelen
\[ \frac{1}{n + 1} \le \frac{\pi^2}{6} - S_n \le \frac{1}{n} \]
kestirimiyle uyumludur (bkz. Analiz 2). Seri yavaş yakınsar: bir milyon terimden sonra bile hata \(10^{-6}\) mertebesindedir. \(\blacksquare\)
6.3 Boolean Maskeler ve np.where
Karşılaştırma işleçleri de evrensel fonksiyondur ve sonuçları True/False değerlerinden oluşan boolean dizilerdir. Böyle bir maskeyle eleman seçmeyi (x[mask]), koşulları & (ve), | (veya), ~ (değil) ile birleştirmeyi ve mask.sum() ile saymayı NumPy Dizileri bölümünde görmüştük. Maskenin sayaç olarak birkaç kullanımı daha vardır:
import numpy as np
x = np.arange(-3, 4)
mask = x > 0
print(mask)
print(mask.sum(), mask.mean()) # True sayısı ve oranı
print(np.any(x > 2), np.all(x > -5))Çıktı:
[False False False False True True True]
3 0.42857142857142855
True TrueTrue değeri \(1\), False değeri \(0\) sayıldığından mask.sum() koşulu sağlayan eleman sayısını, mask.mean() ise bunların oranını verir: yedi elemanın üçü pozitiftir ve oran \(3/7 \approx 0{,}43\)’tür. np.any en az bir elemanın, np.all bütün elemanların koşulu sağlayıp sağlamadığını söyler.
NumPy Dizileri bölümünde and, or ve not sözcüklerinin dizilerle çalışmadığını görmüştük. Parantezi unutmak da aynı hataya götürür. & işlecinin önceliği karşılaştırmalardan yüksek olduğu için x > -2 & x < 2 ifadesi x > (-2 & x) < 2 diye okunur. Bu bir zincirleme karşılaştırmadır ve Python onu içeride and ile çözer. Yani yine bir dizinin bütün olarak doğru olup olmadığını sormuş oluruz:
import numpy as np
x = np.arange(-3, 4)
try:
x[x > -2 & x < 2] # parantez yok: -2 & x önce hesaplanır
except ValueError as err:
print("Hata:", str(err).split(".")[0]) # mesajın ilk cümlesi
print(x[(x > -2) & (x < 2)]) # doğrusu: her koşul parantez içindeÇıktı:
Hata: The truth value of an array with more than one element is ambiguous
[-1 0 1]Maskeyle eleman seçmek diziyi kısaltır. Çoğu zaman ise dizinin şekli korunsun, koşulu sağlayan ve sağlamayan elemanlar farklı değerler alsın isteriz. Parçalı tanımlı fonksiyonlar tam olarak böyledir.
Tanım 6.3 (np.where ile Eleman Eleman Seçim) \(c\) bir boolean dizisi, \(a\) ve \(b\) aynı şekle sahip diziler ya da sayılar olsun. np.where(c, a, b) çağrısı, her \(i\) için
\[ z_i = \begin{cases} a_i, & c_i \text{ doğru ise}, \\ b_i, & \text{diğer durumlarda} \end{cases} \]
elemanlarından oluşan \(z\) dizisini verir. Tek argümanla çağrılan np.where(c) ise \(c\)’nin doğru olduğu indeksleri döndürür.
Yani np.where, “eğer … ise … değilse …” kalıbının eleman eleman çalışan biçimidir. Örneğin
\[ f(x) = \begin{cases} 1 - x^2, & |x| \le 1, \\ 0, & \text{diğer durumlarda} \end{cases} \]
fonksiyonu tek satırda hesaplanır:
import numpy as np
x = np.linspace(-2, 2, 9)
f = np.where(np.abs(x) <= 1, 1 - x**2, 0.0)
print(x)
print(f)
print(np.where(x > 0)) # tek argüman: koşulun doğru olduğu indeksler
print(np.where(x > 0, x, -x) == np.abs(x))Çıktı:
[-2. -1.5 -1. -0.5 0. 0.5 1. 1.5 2. ]
[0. 0. 0. 0.75 1. 0.75 0. 0. 0. ]
(array([5, 6, 7, 8]),)
[ True True True True True True True True True]Son satır, np.where(x > 0, x, -x) ifadesinin tam olarak mutlak değer olduğunu doğruluyor. Tek argümanlı çağrı bir demet döndürür. Dizinin her ekseni için bir indeks dizisi vardır, burada eksen tek olduğu için demetin tek elemanı var.
Örnek 6.4 (Parçalı Fonksiyonun İntegrali) Yukarıdaki \(f\) fonksiyonunun \([-2, 2]\) aralığındaki integralini, aralığı \(n = 4000\) eşit parçaya bölen orta nokta toplamıyla hesaplayın ve kesin değer \(4/3\) ile karşılaştırın.
Çözüm
\(h = (b - a)/n\) olmak üzere alt aralıkların orta noktaları \(m_i = a + (i + 1/2)h\), \(i = 0, 1, \ldots, n - 1\) noktalarıdır. Orta nokta toplamı \(h \sum_i f(m_i)\)’dir (Riemann toplamı için bkz. Analiz 2). Orta noktaların hepsi bir dizide, fonksiyon değerleri np.where ile, toplam da sum ile bulunur:
import numpy as np
n = 4000
a, b = -2.0, 2.0
h = (b - a) / n
m = a + h * (np.arange(n) + 0.5) # alt aralıkların orta noktaları
f = np.where(np.abs(m) <= 1, 1 - m**2, 0.0)
approx = h * f.sum()
print(approx, 4 / 3, abs(approx - 4 / 3))Çıktı:
1.3333335 1.3333333333333333 1.66666666689963e-07Kesin değer \(\int_{-1}^{1} (1 - x^2)\, dx = 4/3\)’tür. Hata yaklaşık \(1{,}67 \cdot 10^{-7}\)’dir. Hata yalnız \(f\)’nin sıfır olmadığı \([-1, 1]\) aralığından gelir. Orta nokta kuralının bu aralıktaki hatası, aralığın uzunluğu \(L\) olmak üzere \(L h^2 |f''|/24\) formülüyle verilir. Burada \(L = 2\), \(h^2 = 10^{-6}\) ve \(|f''| = 2\) olduğundan bu formül tam olarak \(1{,}67 \cdot 10^{-7}\) değerini verir. Kırılma noktaları \(\pm 1\) ızgaranın düğümlerine denk geldiği için fazladan bir hata da oluşmadı. \(\blacksquare\)
np.where(c, a, b) çağrılmadan önce a ve b dizileri bütünüyle hesaplanır; seçim ondan sonra yapılır. np.where(x >= 0, np.sqrt(x), 0.0) yazınca karekök negatif sayılar için de alınır. Sonuç doğru çıkar ama bir uyarı basılır. Aşağıda np.errstate(invalid="raise") uyarıyı hataya çevirerek bunu görünür kılıyor. Geçersiz değerler fonksiyona hiç gönderilmemelidir: ya girdi önceden düzeltilir (np.maximum(x, 0)) ya da ufunc’ın where= parametresi kullanılır.
import numpy as np
x = np.array([-4.0, -1.0, 0.0, 1.0, 4.0])
with np.errstate(invalid="raise"): # uyarıyı hataya çevir
try:
y = np.where(x >= 0, np.sqrt(x), 0.0)
except FloatingPointError as err:
print("Hata:", err)
y1 = np.sqrt(np.maximum(x, 0.0))
y2 = np.sqrt(x, out=np.zeros_like(x), where=x >= 0)
print(y1)
print(y2)Çıktı:
Hata: invalid value encountered in sqrt
[0. 0. 0. 1. 2.]
[0. 0. 0. 1. 2.]Maskelerin güzel bir uygulaması, bir olayın olasılığını rastgele denemelerle kestirmektir. \([0, 1) \times [0, 1)\) karesine düzgün dağılımla nokta atarsak bir noktanın \(x^2 + y^2 \le 1\) çeyrek dairesine düşme olasılığı alanların oranı, yani \(\pi/4\) olur. Öyleyse içeri düşen noktaların oranının 4 katı \(\pi\) için bir kestirimdir:
import numpy as np
rng = np.random.default_rng(2024)
n = 1000
P = rng.random((n, 2)) # birim karede n nokta, şekil (n, 2)
inside = (P**2).sum(axis=1) <= 1.0 # çeyrek dairenin içinde mi?
print(inside[:8])
print(inside.sum(), 4 * inside.mean())
for n in [10**4, 10**5, 10**6]:
P = rng.random((n, 2))
print(n, 4 * np.mean((P**2).sum(axis=1) <= 1.0))Çıktı:
[ True True False True True True True True]
776 3.104
10000 3.162
100000 3.1472
1000000 3.141404rng.random((n, 2)) her satırı bir nokta olan \(n \times 2\) bir dizi üretir. (P**2).sum(axis=1) her noktanın \(x^2 + y^2\) değerini verir, karşılaştırma bunu maskeye çevirir ve mean içerideki noktaların oranını hesaplar. Noktalar üzerinde hiçbir döngü yok; tek döngü, denenen \(n\) değerleri üzerindedir. Kestirim \(n\) büyüdükçe \(\pi\)’ye yaklaşır, ama yavaş: hata kabaca \(1/\sqrt{n}\) ile küçülür. Rastgele sayı üretecini, bu hata hızını ve Monte Carlo yöntemlerini Rastgele Sayılar ve Monte Carlo Yöntemleri bölümünde ayrıntılı ele alacağız.
default_rng(2024) ile üretilen ilk 1000 nokta. x² + y² ≤ 1 maskesi 776 noktada doğrudur (mavi). Çeyrek dairenin alanı π/4 olduğundan 4 · 776/1000 = 3,104 değeri π için bir kestirimdir.Örnek 6.5 (Kaydırılmış Maskelerle İkiz Asallar) \(p\) ile \(p + 2\) sayılarının ikisi de asalsa bu çifte ikiz asal çifti denir. İkisi de \(10^6\)’dan küçük olan ikiz asal çiftlerinin sayısını, Eratosthenes kalburunun maskesini kaydırarak döngüsüz bulun.
Çözüm
Kalbur maskesini NumPy Dizileri bölümündeki gibi kurarız; yalnız \(N\) bu kez \(10^6\)’dır. is_prime[p] değeri \(p\)’nin asal olup olmadığını söyler. Python döngüsü yalnız \(p \le 1000\) için döner; asıl iş, yani yaklaşık 2,1 milyon silme, vektörel dilim atamalarında yapılır.
Yeni olan kısım ikiz asallardır. is_prime[:-2] dilimi \(p = 0, 1, \ldots, N - 3\) için \(p\)’nin, is_prime[2:] dilimi ise aynı konumda \(p + 2\)’nin asal olup olmadığını söyler. İki dilim aynı uzunluktadır, bu yüzden & ile birleşimleri tam olarak \(p\) ile \(p + 2\)’nin ikisinin de asal olduğu konumlarda True olur. np.flatnonzero, tek argümanlı np.where gibi, maskenin True olduğu indeksleri verir.
import numpy as np
N = 10**6
is_prime = np.ones(N, dtype=bool) # 0, 1, ..., N - 1 için bayrak
is_prime[:2] = False # 0 ve 1 asal değil
for p in range(2, int(N**0.5) + 1):
if is_prime[p]:
is_prime[p * p::p] = False # p*p, p*p + p, ... siliniyor
print(is_prime.sum()) # N'den küçük asalların sayısı
twin = is_prime[:-2] & is_prime[2:] # twin[p]: p ve p + 2 asal mı?
first = np.flatnonzero(twin) # her çiftin küçük elemanı
print(first[:8], first[-2:])
print(twin.sum())Çıktı:
78498
[ 3 5 11 17 29 41 59 71] [999611 999959]
8169\(10^6\)’dan küçük \(78498\) asal ve \(8169\) ikiz asal çifti vardır. İlk sekiz çift, Sınıflar, Hata Yakalama ve Dosyalar bölümünde 80’e kadar deneme bölmesiyle bulduğumuz \((3, 5), (5, 7), \ldots, (71, 73)\) çiftleridir; sonuncusu \((999959, 999961)\)’dir. Orada her sayı ayrı ayrı sınanıyordu. Burada bir milyon sayının hepsi tek bir & işlemiyle sınandı. İkiz asal çiftlerinin sonsuz sayıda olup olmadığı hâlâ açık bir sorudur. \(\blacksquare\)
6.4 Broadcasting
Şimdiye kadar ya aynı şekildeki iki diziyi ya da bir diziyle bir sayıyı işledik. A + 5 yazınca 5 sayısı A’nın her elemanına eklendi. Şekilleri farklı iki dizi için de benzer bir kural vardır:
import numpy as np
A = np.arange(12).reshape(3, 4)
row = np.array([100, 200, 300, 400]) # şekil (4,)
col = np.array([[1000], [2000], [3000]]) # şekil (3, 1)
print(A + 5)
print(A + row)
print(A + col)
try:
A + np.array([1, 2, 3]) # şekil (3,)
except ValueError as err:
print("Hata:", err)Çıktı:
[[ 5 6 7 8]
[ 9 10 11 12]
[13 14 15 16]]
[[100 201 302 403]
[104 205 306 407]
[108 209 310 411]]
[[1000 1001 1002 1003]
[2004 2005 2006 2007]
[3008 3009 3010 3011]]
Hata: operands could not be broadcast together with shapes (3,4) (3,) (4,) şekilli row dizisi matrisin her satırına, (3, 1) şekilli col sütunu ise her sütununa eklendi. (3,) şekilli bir dizi ise eklenemedi. Hangi şekillerin birlikte işlenebileceğini aşağıdaki kurallar belirler.
Tanım 6.4 (Broadcasting Kuralları) Şekilleri \((a_1, \ldots, a_p)\) ve \((b_1, \ldots, b_q)\) olan iki dizi bir ufunc ile birlikte işlenirken şu adımlar izlenir:
- Kısa olan şeklin soluna 1’ler eklenerek iki dizinin eksen sayısı eşitlenir.
- Şekiller sağdan sola eksen eksen karşılaştırılır. Her eksende iki uzunluk ya eşit olmalı ya da biri 1 olmalıdır. Aksi hâlde işlem
ValueErrorhatası verir. - Sonucun her eksendeki uzunluğu iki uzunluğun büyüğüdür. Uzunluğu 1 olan eksen, o uzunluğa kadar kopyalanmış gibi davranır; bu kopyalar bellekte oluşturulmaz.
Bu kurallara broadcasting (yayma) denir. Bir sayı, şekli \(()\) olan bir dizi gibi işlem görür.
Yani küçük dizi, büyük dizinin şekline “yayılır”. Uzunluğu 1 olan bir eksen istenen uzunluğa kadar uzatılır, eksik eksenler sola eklenir. Yukarıda (3, 4) ile (4,) şekilleri uyumluydu: (4,) önce (1, 4) olur, sonra satır üç kez kopyalanmış gibi davranır. (3, 4) ile (3,) ise uyumsuzdur, çünkü (3,) dizisi (1, 3) olur ve son eksende 4 ile 3 çatışır. İki dizi aynı anda da yayılabilir. (3, 1) ile (4,) şekillerinden (3, 4) şekilli bir sonuç çıkar:
a + b çıktısındaki (3, 4) şekilli dizidir.import numpy as np
a = np.array([[0], [10], [20]]) # şekil (3, 1)
b = np.arange(4) # şekil (4,)
print(a + b)
print(np.broadcast_shapes((3, 1), (4,)))
print(np.broadcast_shapes((5, 3), (3,)))
print(np.broadcast_shapes((8, 1, 6, 1), (7, 1, 5)))Çıktı:
[[ 0 1 2 3]
[10 11 12 13]
[20 21 22 23]]
(3, 4)
(5, 3)
(8, 7, 6, 5)np.broadcast_shapes bir hesap yapmadan yalnız sonucun şeklini söyler. Kuralların adım adım uygulanışı şöyledir:
| Birinci | İkinci | Hizalanmış şekiller | Sonuç |
|---|---|---|---|
| (3, 1) | (4,) | (3, 1) ve (1, 4) | (3, 4) |
| (5, 3) | (3,) | (5, 3) ve (1, 3) | (5, 3) |
| (8, 1, 6, 1) | (7, 1, 5) | (8, 1, 6, 1) ve (1, 7, 1, 5) | (8, 7, 6, 5) |
| (3, 4) | (3,) | (3, 4) ve (1, 3) | hata: 4 ile 3 |
Broadcasting’den yararlanmak için çoğu zaman tek eksenli bir diziyi sütun ya da satır şekline çevirmek gerekir. Bunun için dizinin indeksine np.newaxis (ya da onun kısaltması olan None) yazılır. Bu, o konuma uzunluğu 1 olan yeni bir eksen ekler:
import numpy as np
x = np.array([1, 2, 3])
print(x.shape, x[:, np.newaxis].shape, x[np.newaxis, :].shape)
print(x[:, None].shape, x.reshape(-1, 1).shape)
print(x[:, None] - x) # D[i, j] = x[i] - x[j]Çıktı:
(3,) (3, 1) (1, 3)
(3, 1) (3, 1)
[[ 0 -1 -2]
[ 1 0 -1]
[ 2 1 0]]x[:, None] - x ifadesinde (3, 1) ile (3,) şekilleri (3, 3) şekline yayılır ve \(D_{ij} = x_i - x_j\) farklarının tablosu çıkar. Bu “sütun eksi satır” kalıbı, bütün çiftler üzerinden yapılan hesapların temelidir.
Örnek 6.6 (Çarpım Tablosu) \(1, 2, \ldots, 9\) sayılarının \(9 \times 9\) çarpım tablosunu broadcasting ile kurun. Tablonun elemanları toplamının \((1 + 2 + \cdots + 9)^2 = 2025\) olduğunu ve tablonun simetrik olduğunu doğrulayın.
Çözüm
\(n\) dizisini sütun yaparız: n[:, None] (9, 1), n ise (9,) şeklindedir. Çarpımları (9, 9) şekline yayılır ve \(T_{ij} = n_i n_j\) olur.
import numpy as np
n = np.arange(1, 10) # şekil (9,)
T = n[:, None] * n # (9, 1) ile (9,) -> (9, 9)
print(T[:4, :6])
print(T.shape, T.sum(), n.sum()**2)
print(np.array_equal(T, np.outer(n, n)), np.array_equal(T, T.T))Çıktı:
[[ 1 2 3 4 5 6]
[ 2 4 6 8 10 12]
[ 3 6 9 12 15 18]
[ 4 8 12 16 20 24]]
(9, 9) 2025 2025
True True\[ \sum_{i,j} n_i n_j = \Big(\sum_i n_i\Big)\Big(\sum_j n_j\Big) = 45^2 = 2025 \]
olduğundan iki toplam eşit çıktı. np.outer aynı tabloyu kuran hazır fonksiyondur. Simetri \(T = T^{\mathsf T}\) eşitliğidir; T.T matrisin transpozudur. \(\blacksquare\)
Örnek 6.7 (Bir Veri Matrisini Standartlaştırma) Her satırı bir ölçüm, her sütunu bir değişken olan \(6 \times 3\) bir \(X\) veri matrisinde her sütundan o sütunun ortalamasını çıkarıp sonucu sütunun standart sapmasına bölün. Bu işleme standartlaştırma, sonuçlara \(z\)-değerleri denir. Elde edilen \(Z\) matrisinin her sütununun ortalamasının 0, standart sapmasının 1 olduğunu doğrulayın.
Çözüm
Sütun ortalamaları X.mean(axis=0) ile gelir ve (3,) şeklindedir. (6, 3) şekilli \(X\)’ten çıkarılınca (1, 3) şekline, oradan da (6, 3) şekline yayılır. Yani her satırdan aynı ortalama vektörü çıkarılır. Standart sapmalara bölme de aynı biçimde çalışır:
\[ Z_{ij} = \frac{X_{ij} - \mu_j}{\sigma_j}. \]
Örnek veriyi rng.normal ile üretiyoruz. Bu çağrı normal dağılımlı rastgele sayılar verir: loc ortalamayı, scale yayılımı (standart sapmayı) belirler. Burada ikisi de (3,) şekilli birer dizidir ve size=(6, 3) şekline yayılır. Böylece her sütun kendi ortalaması ve yayılımıyla üretilir.
import numpy as np
rng = np.random.default_rng(7)
X = rng.normal(loc=[10.0, 50.0, -3.0], scale=[1.0, 8.0, 0.5],
size=(6, 3)) # 6 ölçüm, 3 değişken
mu = X.mean(axis=0) # şekil (3,): sütun ortalamaları
sigma = X.std(axis=0) # şekil (3,): sütun std sapmaları
Z = (X - mu) / sigma # (6, 3) ile (3,) -> (6, 3)
print(np.round(X, 2))
print(np.round(Z, 3))
print(np.allclose(Z.mean(axis=0), 0), np.allclose(Z.std(axis=0), 1))Çıktı:
[[10. 52.39 -3.14]
[ 9.11 46.36 -3.5 ]
[10.06 60.72 -3.25]
[ 9.38 53.92 -2.82]
[10.11 42.56 -3.01]
[10.7 39.25 -3.23]]
[[ 0.211 0.439 0.097]
[-1.512 -0.391 -1.623]
[ 0.325 1.587 -0.426]
[-0.99 0.65 1.61 ]
[ 0.413 -0.915 0.684]
[ 1.553 -1.37 -0.343]]
True TrueSütunların ortalamaları (yaklaşık 10, 50 ve \(-3\)) ve yayılımları çok farklıydı. Standartlaştırmadan sonra üçü de aynı ölçeğe geldi. Ortalamalar yuvarlama hatası kadar sıfıra yakın olduğundan karşılaştırmayı == ile değil np.allclose ile yaptık. \(\blacksquare\)
Broadcasting hata vermesi gereken bazı yerlerde sessizce çalışır. En sık görülen durum, (n,) şekilli bir vektörle (n, 1) şekilli bir sütunun farkıdır. İkisi de “n elemanlı bir vektör” gibi görünür, ama farkları (n, n) şekilli bir tablo olur:
import numpy as np
y = np.array([1.0, 2.0, 3.0, 4.0]) # şekil (4,)
y_hat = np.array([[1.1], [1.9], [3.2], [3.8]]) # şekil (4, 1)
r = y - y_hat
print(r.shape) # (4,) beklenirken (4, 4)
r = y - y_hat.ravel() # ya da y_hat[:, 0]
print(r.shape, r)Çıktı:
(4, 4)
(4,) [-0.1 0.1 -0.2 0.2]Bir hesabın sonucu beklenmedik biçimde büyük çıkıyorsa ilk iş ara sonuçların .shape değerlerini yazdırmaktır.
İki değişkenli bir fonksiyonu bir ızgarada hesaplamak broadcasting’in en doğal kullanımıdır. \(x\) değerleri satır, \(y\) değerleri sütun yapılırsa \(F_{ij} = f(x_j, y_i)\) tablosu tek ifadeyle çıkar. Aynı tabloyu np.meshgrid ile de kurabiliriz. Bu fonksiyon iki koordinat dizisini tam ızgara şekline kopyalar; Matplotlib ile Grafik Çizimi bölümündeki eş yükselti eğrileri ve yüzeyler bu dizileri kullanır.
import numpy as np
x = np.linspace(-1, 1, 5) # şekil (5,)
y = np.linspace(0, 1, 3) # şekil (3,)
F = x[np.newaxis, :]**2 + y[:, np.newaxis] # şekil (3, 5)
print(F)
X, Y = np.meshgrid(x, y) # ikisi de şekil (3, 5)
print(X.shape, np.array_equal(F, X**2 + Y))Çıktı:
[[1. 0.25 0. 0.25 1. ]
[1.5 0.75 0.5 0.75 1.5 ]
[2. 1.25 1. 1.25 2. ]]
(3, 5) TrueBurada \(f(x, y) = x^2 + y\)’dir. Tablonun satırları \(y\)’ye, sütunları \(x\)’e karşılık gelir; meshgrid de varsayılan olarak bu düzeni kullanır. Broadcasting sürümü (3, 5) şekilli iki tam diziyi bellekte oluşturmaz, bu yüzden büyük ızgaralarda daha az bellek harcar.
Örnek 6.8 (Bir Çemberin İçindeki Kafes Noktaları) \(x^2 + y^2 \le r^2\) eşitsizliğini sağlayan \((x, y)\) tamsayı çiftlerinin sayısı \(N(r)\) olsun. \(r = 10, 100, 1000\) için \(N(r)\)’yi broadcasting ile hesaplayın ve \(\pi r^2\) ile karşılaştırın.
Çözüm
Kontrol Yapıları ve Fonksiyonlar bölümünde \(N(r)\)’yi iç içe iki döngüyle \(r = 20\)’ye kadar saymıştık. Şimdi döngülerin yerini bir tablo alır. Hem \(x\) hem \(y\), \(-r, \ldots, r\) tamsayılarından birini alır. \(k\) bu tamsayıların dizisi olsun. k[:, None]**2 + k**2 ifadesi \((2r + 1) \times (2r + 1)\) şekilli bir tabloda bütün \(x^2 + y^2\) değerlerini verir. Bunu \(r^2\) ile karşılaştırıp maskedeki True değerlerini sayarız:
import numpy as np
for r in [10, 100, 1000]:
k = np.arange(-r, r + 1) # şekil (2r + 1,)
inside = k[:, None]**2 + k**2 <= r**2 # şekil (2r+1, 2r+1)
N = inside.sum()
print(f"r = {r:>4} N(r) = {N:>7} pi r^2 = {np.pi * r**2:12.2f}"
f" fark = {N - np.pi * r**2:8.2f}")Çıktı:
r = 10 N(r) = 317 pi r^2 = 314.16 fark = 2.84
r = 100 N(r) = 31417 pi r^2 = 31415.93 fark = 1.07
r = 1000 N(r) = 3141549 pi r^2 = 3141592.65 fark = -43.65Orada birim karelerle gerekçelendirdiğimiz \(N(r) \approx \pi r^2\) yaklaşımı çok daha büyük \(r\) değerlerinde de tutuyor: fark, \(\pi r^2\)’nin yanında çok küçüktür. Farkın büyüme hızı Gauss çember problemi olarak bilinen ünlü bir açık sorudur. \(r = 1000\) için tablo yaklaşık 4 milyon elemanlıdır ve yine de bir anda hesaplanır. \(\blacksquare\)
6.5 Uygulama: Uzaklık Matrisi
Broadcasting’in en çok kullanılan uygulamalarından biri, bir nokta kümesindeki bütün çiftlerin uzaklıklarını tek seferde hesaplamaktır. Düzlemde \(m\) nokta \(m \times 2\) şekilli bir \(P\) dizisinde, her satır bir nokta olacak biçimde tutulsun. Öklid uzaklığı (bkz. Analiz 4)
\[ D_{ij} = \lVert P_i - P_j \rVert = \sqrt{(P_{i1} - P_{j1})^2 + (P_{i2} - P_{j2})^2} \]
formülüyle verilir. Burada \(P_{i1}\) ve \(P_{i2}\), \(i\)’inci noktanın iki koordinatıdır. P[:, None, :] (m, 1, 2), P[None, :, :] ise (1, m, 2) şeklindedir. Farkları (m, m, 2) şekline yayılır: diff[i, j] dizisi \(P_i - P_j\) vektörüdür. Son eksen üzerinden kareler toplanıp karekök alınınca (m, m) şekilli uzaklık matrisi çıkar.
import numpy as np
from scipy.spatial.distance import cdist
rng = np.random.default_rng(5)
P = np.round(10 * rng.random((6, 2)), 1) # 6 nokta, şekil (6, 2)
diff = P[:, None, :] - P[None, :, :] # şekil (6, 6, 2)
D = np.sqrt((diff**2).sum(axis=-1)) # şekil (6, 6)
print(P)
print(np.round(D, 2))
print(np.allclose(D, D.T), np.allclose(D, cdist(P, P)))
E = D.copy()
np.fill_diagonal(E, np.inf) # kendisi aday olmasın
print(E.argmin(axis=1)) # en yakın komşularÇıktı:
[[ 8.1 8.1]
[ 5.2 2.9]
[ 0.5 3.8]
[ 4.1 0.5]
[ 0.5 10. ]
[ 6.5 2.3]]
[[ 0. 5.95 8.73 8.59 7.83 6.02]
[ 5.95 0. 4.79 2.64 8.51 1.43]
[ 8.73 4.79 0. 4.88 6.2 6.18]
[ 8.59 2.64 4.88 0. 10.16 3. ]
[ 7.83 8.51 6.2 10.16 0. 9.76]
[ 6.02 1.43 6.18 3. 9.76 0. ]]
True True
[1 5 1 1 2 1]Matris simetriktir ve köşegeni sıfırdır. SciPy’ın cdist fonksiyonu da aynı sonucu verdi. En yakın komşuyu bulmak için her satırın köşegen dışındaki en küçük elemanının yerini ararız. Köşegeni np.inf yaparsak bir nokta kendisini seçemez ve argmin(axis=1) her nokta için en yakın komşunun indeksini verir.
default_rng(5) ile üretilen altı nokta; her ok bir noktadan en yakın komşusuna gider (1 ile 5 birbirinin en yakın komşusudur). Sağda broadcasting ile hesaplanan uzaklık matrisi D: uzaklık küçüldükçe hücre koyulaşır, her satırda köşegen dışındaki en küçük eleman çerçevelidir. Çerçeveli sütunlar E.argmin(axis=1) sonucu olan [1 5 1 1 2 1] dizisidir.Bu yöntemin bir bedeli vardır: ara dizi diff \(m^2 d\) eleman tutar (\(d\) uzayın boyutudur). \(m = 10^4\) nokta ve \(d = 3\) için bu \(3 \cdot 10^8\) sayı, yani 2,4 GB bellek eder. Vektörizasyon hızı çoğu zaman bellekle satın alır. Çok büyük nokta kümelerinde hesap parçalara bölünür ya da cdist gibi ara dizi kurmayan hazır fonksiyonlar kullanılır.
6.6 Döngü mü, Vektör mü?
Vektörel kodun daha hızlı olduğunu söyledik; şimdi bunu ölçelim. Python’un timeit modülü bir ifadeyi birçok kez çalıştırıp süresini ölçer. Bilgisayarda aynı anda çalışan başka işler ölçümü yalnız yavaşlatabileceği için tekrarların en küçüğünü almak âdettir. Bir milyon sayının karelerinin toplamını üç yolla hesaplayalım:
import timeit
import numpy as np
n = 10**6
x = np.random.default_rng(0).random(n)
xs = x.tolist() # aynı sayılar, Python listesi olarak
def sum_squares_index(a):
total = 0.0
for i in range(len(a)): # dizinin elemanlarına tek tek eriş
total += a[i] * a[i]
return total
def sum_squares_loop(values):
total = 0.0
for v in values:
total += v * v
return total
print(sum_squares_index(x), sum_squares_loop(xs), np.sum(x * x))
def best_time(f, number):
return min(timeit.repeat(f, number=number, repeat=7)) / number
t_index = best_time(lambda: sum_squares_index(x), 1)
t_loop = best_time(lambda: sum_squares_loop(xs), 1)
t_vec = best_time(lambda: np.sum(x * x), 10)
print(f"dizi üzerinde döngü : {t_index * 1e3:6.0f} ms")
print(f"liste üzerinde döngü : {t_loop * 1e3:6.0f} ms")
print(f"vektörel : {t_vec * 1e3:6.1f} ms")Çıktı:
333560.6157428297 333560.6157428297 333560.61574281874
dizi üzerinde döngü : 128 ms
liste üzerinde döngü : 19 ms
vektörel : 1.4 msÜç yöntem de aynı toplamı verdi. Döngünün sonucuyla np.sum’ınki son basamaklarda farklıdır. np.sum sayıları gruplara ayırıp önce grup toplamlarını alır (ikili toplama), bu yüzden yuvarlama hatası daha az birikir. Süreler bu bölümün yazıldığı bilgisayarda ölçülmüştür; sizin makinenizde farklı çıkacaktır, ama oranların mertebesi genellikle benzer kalır. Vektörel toplam liste üzerindeki döngüden yaklaşık 14 kat, dizi üzerindeki döngüden yaklaşık 90 kat hızlıdır.
Dizi üzerindeki döngünün listeden bile yavaş olması şaşırtıcı görünebilir. Bunun nedeni, a[i] ile her erişimde NumPy’ın makine sayısını bir Python nesnesine sarmak zorunda olmasıdır. Döngü her elemanda bu paketlemeyi, tür denetimini ve yorumlayıcının işini tekrarlar. Vektörel kodda ise bu masraf dizi başına bir kez ödenir ve elemanlar üzerindeki döngü derlenmiş kodda, bellekte art arda duran sayılar üzerinde döner. Farkın \(n\) ile nasıl değiştiğini görmek için ölçümü farklı boyutlarda tekrarlayalım. Tabloya bir de np.vectorize ile “vektörleştirilmiş” bir Python fonksiyonu ekleyelim:
import timeit
import numpy as np
def sum_squares_loop(values):
total = 0.0
for v in values:
total += v * v
return total
square = np.vectorize(lambda v: v * v)
rng = np.random.default_rng(0)
print(" n döngü (s) vectorize (s) vektörel (s)")
for n in [10, 100, 10**3, 10**4, 10**5, 10**6]:
x = rng.random(n)
xs = x.tolist()
reps = max(1, 10**5 // n)
times = []
for f in [lambda: sum_squares_loop(xs),
lambda: np.sum(square(x)),
lambda: np.sum(x * x)]:
t = min(timeit.repeat(f, number=reps, repeat=5)) / reps
times.append(t)
print(f"{n:>7} {times[0]:9.1e} {times[1]:11.1e}"
f" {times[2]:11.1e}")Çıktı:
n döngü (s) vectorize (s) vektörel (s)
10 2.6e-07 5.7e-06 1.4e-06
100 2.0e-06 1.2e-05 1.4e-06
1000 2.0e-05 7.6e-05 2.0e-06
10000 2.1e-04 7.0e-04 5.0e-06
100000 2.0e-03 8.0e-03 1.6e-04
1000000 1.9e-02 8.4e-02 1.5e-03np.vectorize içeride yine Python döngüsü çalıştırdığı için düz döngüden de yavaştır. Sayılar makineye göre değişir; eğrilerin biçimi değişmez.Tablo üç şey söylüyor. Birincisi, çok küçük dizilerde (10 eleman) vektörel kodun yaklaşık 1 mikrosaniyelik sabit başlangıç maliyeti döngüden pahalıdır; vektörizasyon büyük dizilerde kazandırır. İkincisi, \(n \ge 10^3\) için vektörel kod yaklaşık 10 ile 40 kat hızlıdır. Vektörel eğrinin \(10^4\) ile \(10^5\) arasındaki sıçraması, x * x ara dizisi için bellek ayırma ve önbellek etkilerinden gelir; ölçümün makineden makineye, hatta çalıştırmadan çalıştırmaya en çok değişen kısmı budur. Üçüncüsü, np.vectorize hiç hızlandırmaz.
np.vectorize(f), tek sayı alan bir Python fonksiyonunu dizilere uygulanabilir hâle getirir. Ama bunu içeride bir Python döngüsüyle yapar. Kodu kısaltır, ama yukarıdaki tabloda görüldüğü gibi düz döngüden bile yavaştır. Gerçek vektörizasyon, hesabı ufunc’lar, dilimler, maskeler ve broadcasting ile yeniden yazmaktır.
6.7 Uygulama: Sayısal Türev
Türev tanımındaki fark bölümü, kaydırılmış dizilerin farkı olarak tek satırda yazılabilir. \([0, 2\pi]\) aralığını \(16\) eşit parçaya bölen \(x_0, x_1, \ldots, x_{16}\) noktalarında \(y_i = \sin x_i\) değerleri bilinsin ve \(h = 2\pi/16\) olsun. np.diff(y) dizisi \(y_{i+1} - y_i\) farklarını verir; bu farklar 17 değil 16 tanedir. y[2:] - y[:-2] ise \(y_{i+1} - y_{i-1}\) farklarıdır; bunlar iç noktalara aittir ve sayıları 15’tir.
import numpy as np
n = 17
x = np.linspace(0, 2 * np.pi, n) # 16 eşit alt aralık
h = x[1] - x[0]
y = np.sin(x)
d_fwd = np.diff(y) / h # n - 1 değer
d_cen = (y[2:] - y[:-2]) / (2 * h) # n - 2 değer, x[1:-1] için
mid = (x[:-1] + x[1:]) / 2 # alt aralıkların orta noktaları
print(round(h, 4), len(d_fwd), len(d_cen))
print(np.abs(d_fwd - np.cos(x[:-1])).max()) # sol uca yazılırsa
print(np.abs(d_fwd - np.cos(mid)).max()) # orta noktaya yazılırsa
print(np.abs(d_cen - np.cos(x[1:-1])).max()) # merkezi farkÇıktı:
0.3927 16 15
0.19383917874071455
0.006289921998797743
0.025504641595567312İleri fark bölümü \((y_{i+1} - y_i)/h\), \(x_i\) noktasındaki türeve ancak \(O(h)\) doğrulukla yaklaşır (bkz. Nümerik Analiz). Aynı sayı \(x_i\) ile \(x_{i+1}\) arasındaki orta noktadaki türeve ise \(O(h^2)\) doğrulukla yaklaşır. Yani hangi değerin hangi noktaya ait olduğu, hesabın kendisi kadar önemlidir. Aynı 16 sayı sol uçlara yazılınca en büyük hata \(0{,}19\), orta noktalara yazılınca \(0{,}0063\) çıktı. Merkezi fark formülü \((y_{i+1} - y_{i-1})/(2h)\) iç noktalarda yine \(O(h^2)\) doğruluk verir (bkz. Nümerik Analiz).
np.diff(y) / h ile bulunan 16 fark bölümü (h ≈ 0,39). Değerler alt aralıkların sol uçlarına yazılınca (boş daireler) cos x eğrisinin yarım adım soluna kayar, en büyük hata 0,19 olur. Aynı sayılar orta noktalara yazılınca (dolu daireler) eğrinin üstüne oturur, en büyük hata 0,0063'e iner.- Noktaları bir
xdizisinde kur vehadımını belirle. Fonksiyon değerlerini tek ifadeyle hesapla:f(x),f(x + h), … ya da ızgaraday = f(x). - Formülü kaydırılmış dizilerle yaz: ızgarada
y[1:],y[:-1],y[2:]gibi dilimlerle, fonksiyon biliniyorsaf(x + h) - f(x - h)gibi ifadelerle. Birden çokhdeğerini aynı anda denemek içinh’yi sütun yap. - Sonucun hangi noktalara ait olduğunu belirle (
x[1:-1], orta noktalar, …). Kesin türev biliniyorsa farkın mutlak değerinin en büyüğünü uygunaxisüzerinden al.
Örnek 6.9 (Merkezi Farkın Hatası ve Adım Boyu) \(f(x) = \sin x\) için merkezi fark formülü \(\big(f(x + h) - f(x - h)\big)/(2h)\) ile \([0, 2\pi]\) aralığındaki 201 noktada türevi \(h = 10^{-1}, 10^{-2}, 10^{-3}\) adımlarıyla hesaplayın. Her \(h\) için en büyük hatayı bulun ve hatanın \(h^2\) ile orantılı küçüldüğünü gösterin.
Çözüm
1. adım. x (201,) şekilli ızgaradır. h ise (3, 1) şeklinde bir sütun olarak kurulur.
2. adım. np.sin(x + h) ifadesinde (3, 1) ile (201,) şekilleri (3, 201) şekline yayılır. Böylece her satır bir \(h\) değeri için 201 noktadaki fark bölümlerini tutar.
3. adım. Sonuç x noktalarının kendisine aittir. Kesin türev \(\cos x\)’ten farkın mutlak değeri axis=1 boyunca, yani her \(h\) için ayrı ayrı en büyüğe indirgenir.
import numpy as np
x = np.linspace(0, 2 * np.pi, 201) # şekil (201,)
h = np.array([1e-1, 1e-2, 1e-3])[:, None] # şekil (3, 1)
D = (np.sin(x + h) - np.sin(x - h)) / (2 * h) # şekil (3, 201)
err = np.abs(D - np.cos(x)).max(axis=1) # şekil (3,)
for hk, ek in zip(h.ravel(), err):
print(f"h = {hk:.0e} en büyük hata = {ek:.3e}"
f" h^2/6 = {hk**2 / 6:.3e}")
print(err[:-1] / err[1:])Çıktı:
h = 1e-01 en büyük hata = 1.666e-03 h^2/6 = 1.667e-03
h = 1e-02 en büyük hata = 1.667e-05 h^2/6 = 1.667e-05
h = 1e-03 en büyük hata = 1.667e-07 h^2/6 = 1.667e-07
[99.95051153 99.99943903]Merkezi farkın hatası \(\frac{h^2}{6} f'''(\xi)\)’dir ve \(|f'''| = |\cos| \le 1\) olduğundan en büyük hata \(h^2/6\) civarında beklenir. Tablo bunu çok iyi doğruluyor. \(h\) on kat küçülünce hata yüz kat küçülüyor (son satırdaki oranlar), yani yöntem ikinci mertebedendir. \(h\)’yi çok küçültmenin bir sınırı da vardır, çünkü bir noktadan sonra yuvarlama hatası baskın çıkar. Bunu Sayısal Türev ve İntegral bölümünde inceleyeceğiz. \(\blacksquare\)
6.8 Alıştırmalar
Aşağıdaki alıştırmalarda döngü yerine ufunc’lar, indirgemeler, maskeler ve broadcasting kullanacağız. Her çözüm, kodu ve çalıştırılmış çıktısıyla birlikte verilmiştir; önce kendiniz deneyip sonra karşılaştırın.
Alıştırma 6.1 (Wallis Çarpımı) Wallis çarpımı \(\displaystyle \frac{\pi}{2} = \prod_{k=1}^{\infty} \frac{4k^2}{4k^2 - 1}\) eşitliğini söyler (bkz. Analiz 2). \(W_n = 2\prod_{k=1}^{n} \frac{4k^2}{4k^2 - 1}\) kısmi çarpımlarını \(n = 10^5\)’e kadar np.cumprod ile hesaplayın ve \(\pi - W_n\) farkının \(n\) ile nasıl küçüldüğünü bulun.
Çözüm
k bir tamsayı dizisidir, ama / bölmesi tamsayı dizilerinde de float sonuç verir; bu yüzden çarpanlar doğru hesaplanır. np.cumprod bütün kısmi çarpımları verir; \(W_n\) değeri W[n - 1] konumundadır.
import numpy as np
k = np.arange(1, 10**5 + 1) # 1, 2, ..., 10^5
factors = 4 * k**2 / (4 * k**2 - 1)
W = 2 * np.cumprod(factors) # W[n - 1]: ilk n çarpan
for n in [10, 100, 1000, 10**4, 10**5]:
gap = np.pi - W[n - 1]
print(f"n = {n:>6} W_n = {W[n - 1]:.8f} "
f"fark = {gap:.3e} n * fark = {n * gap:.4f}")Çıktı:
n = 10 W_n = 3.06770381 fark = 7.389e-02 n * fark = 0.7389
n = 100 W_n = 3.13378749 fark = 7.805e-03 n * fark = 0.7805
n = 1000 W_n = 3.14080775 fark = 7.849e-04 n * fark = 0.7849
n = 10000 W_n = 3.14151412 fark = 7.853e-05 n * fark = 0.7853
n = 100000 W_n = 3.14158480 fark = 7.854e-06 n * fark = 0.7854Son sütun \(n \cdot (\pi - W_n)\) değerinin \(0{,}7854 \approx \pi/4\) sayısına yaklaştığını gösteriyor. Yani \(\pi - W_n \approx \dfrac{\pi}{4n}\)’dir. Wallis çarpımı da Basel serisi gibi \(1/n\) hızıyla, yavaş yakınsar: \(10^5\) çarpandan sonra bile hata yaklaşık \(8 \cdot 10^{-6}\)’dır. \(\blacksquare\)
Alıştırma 6.2 (Sütunların En Büyük Elemanları) rng = np.random.default_rng(3) üreteciyle rng.integers(0, 100, size=(5, 4)) matrisini kurun. Her sütunun en büyük elemanını ve bu elemanın bulunduğu satırın indeksini döngü kullanmadan bulun.
Çözüm
Sütun boyunca indirgeme axis=0 ile yapılır. A.max(axis=0) en büyük değerleri, A.argmax(axis=0) bunların satır indekslerini verir. Bir sütunda en büyük değer birden çok kez geçerse argmax ilkinin indeksini döndürür. Son satırda fancy indeksleme yaparak (bkz. NumPy Dizileri) \(A_{i_j j}\) elemanlarını alıyor ve sonucun en büyük değerlerle aynı olduğunu doğruluyoruz.
import numpy as np
rng = np.random.default_rng(3)
A = rng.integers(0, 100, size=(5, 4))
print(A)
print(A.max(axis=0)) # her sütunun en büyük elemanı
print(A.argmax(axis=0)) # bu elemanların satır indeksleri
cols = np.arange(A.shape[1])
print(A[A.argmax(axis=0), cols])Çıktı:
[[81 8 17 23]
[18 80 86 58]
[ 3 9 33 43]
[62 47 26 15]
[69 73 3 11]]
[81 80 86 58]
[0 1 1 1]
[81 80 86 58]Sütunların en büyükleri \(81, 80, 86, 58\)’dir. İlki 0 numaralı satırda, diğer üçü 1 numaralı satırdadır. \(\blacksquare\)
Alıştırma 6.3 (Satırları Normalleştirme) Elemanları negatif olmayan
\[ W = \begin{pmatrix} 2 & 1 & 1 \\ 0 & 3 & 1 \\ 1 & 1 & 6 \end{pmatrix} \]
matrisinin her satırını kendi toplamına bölerek satır toplamları 1 olan bir \(P\) matrisi kurun. Böyle matrislere satır stokastik matris denir; Markov zincirlerinin geçiş matrisleri bu türdendir.
Çözüm
Satır toplamları W.sum(axis=1) ile bulunur, ama bu (3,) şekilli bir dizidir. (3, 3) şekilli \(W\)’yi buna bölersek (3,) dizisi (1, 3) şekline yayılır ve \(j\)’inci sütun \(s_j\)’ye bölünür; istediğimiz bu değildir. Toplamı keepdims=True ile (3, 1) şeklinde bir sütun olarak tutarsak \(i\)’inci satır \(s_i\)’ye bölünür:
import numpy as np
W = np.array([[2.0, 1.0, 1.0],
[0.0, 3.0, 1.0],
[1.0, 1.0, 6.0]])
s = W.sum(axis=1, keepdims=True) # şekil (3, 1)
P = W / s # her satır kendi toplamına bölünür
print(s.ravel())
print(P)
print(P.sum(axis=1))
print(W / W.sum(axis=1)) # yanlış: j. sütun s_j'ye bölünürÇıktı:
[4. 4. 8.]
[[0.5 0.25 0.25 ]
[0. 0.75 0.25 ]
[0.125 0.125 0.75 ]]
[1. 1. 1.]
[[0.5 0.25 0.125]
[0. 0.75 0.125]
[0.25 0.25 0.75 ]]Satır toplamları \(4, 4, 8\)’dir. \(P\)’nin her satırının toplamı 1 çıktı. Son çıktı yanlış yolun sonucudur. Hata vermedi, çünkü matris karedir ve şekiller uyumludur. Ama sonuç, örneğin ilk satırı \((0{,}5;\ 0{,}25;\ 0{,}125)\), istediğimiz matris değildir. \(\blacksquare\)
Alıştırma 6.4 (Sonucun Şeklini Önceden Bulmak) A = np.ones((6, 1, 4)) ve B = np.arange(5).reshape(5, 1) dizileri için A + B toplamının şeklini broadcasting kurallarıyla elle bulun, sonra kodla doğrulayın.
Çözüm
1. kural. B’nin şekli (5, 1) iki eksenlidir; soluna 1 eklenince (1, 5, 1) olur.
2. kural. Sağdan sola karşılaştırırız: son eksende 4 ile 1, ortada 1 ile 5, ilk eksende 6 ile 1. Her çiftte ya uzunluklar eşit ya da biri 1’dir, yani şekiller uyumludur.
3. kural. Her eksende büyük uzunluk alınır: sonuç (6, 5, 4) şeklindedir.
import numpy as np
A = np.ones((6, 1, 4))
B = np.arange(5).reshape(5, 1)
C = A + B
print(C.shape, np.broadcast_shapes(A.shape, B.shape))
print(C[0]) # C[i, j, k] = 1 + B[j, 0]Çıktı:
(6, 5, 4) (6, 5, 4)
[[1. 1. 1. 1.]
[2. 2. 2. 2.]
[3. 3. 3. 3.]
[4. 4. 4. 4.]
[5. 5. 5. 5.]]C[i, j, k] = A[i, 0, k] + B[j, 0] = 1 + j olduğundan her C[i] dilimi, \(j\)’inci satırı \(1 + j\) olan \(5 \times 4\) bir matristir. \(\blacksquare\)
Alıştırma 6.5 (İşaret Değişimiyle Kök Aralığı) \(f(x) = x^3 - 2x - 5\) fonksiyonunu \([-3, 3]\) aralığında \(0{,}25\) adımlı ızgarada hesaplayın ve işaretin değiştiği alt aralıkları döngü kullanmadan bulun.
Çözüm
np.sign(f) her değerin işaretini (\(-1\), \(0\) ya da \(1\)) verir. Komşu iki değerin işaretleri farklıysa çarpımları negatiftir: s[:-1] * s[1:] < 0 maskesi \([x_i, x_{i+1}]\) aralığında işaret değişimi olan \(i\) indekslerinde doğrudur. np.flatnonzero bu indeksleri listeler.
import numpy as np
x = np.linspace(-3, 3, 25) # adım 0.25
f = x**3 - 2 * x - 5
s = np.sign(f)
idx = np.flatnonzero(s[:-1] * s[1:] < 0)
print(idx)
for i in idx:
print(x[i], x[i + 1], f[i], f[i + 1])Çıktı:
[20]
2.0 2.25 -1.0 1.890625İşaret yalnız \([2;\ 2{,}25]\) aralığında değişir: \(f(2) = -1 < 0\), \(f(2{,}25) = 1{,}890625 > 0\). \(f\) sürekli olduğundan Ara Değer Teoremi (bkz. Nümerik Analiz) gereği bu aralıkta bir kök vardır (bu kök yaklaşık \(2{,}0946\)’dır). Böyle bir aralık, ikiye bölme gibi kök bulma yöntemleri için iyi bir başlangıçtır; bkz. Kök Bulma ve Optimizasyon. \(\blacksquare\)
Alıştırma 6.6 (Izgarada Yerel Maksimumlar) \(f(x) = e^{-x/5} \sin x\) fonksiyonunu \([0, 20]\) aralığında \(0{,}01\) adımlı ızgarada hesaplayın. İki komşusundan da büyük olan iç ızgara noktalarını döngü kullanmadan bulun ve \(f'(x) = 0\) koşulundan gelen kesin yerel maksimum noktalarıyla karşılaştırın.
Çözüm
1. adım. Önce kesin yerleri bulalım. \(f'(x) = e^{-x/5}\big(\cos x - \tfrac{1}{5}\sin x\big)\) olduğundan kritik noktalar \(\tan x = 5\), yani \(x = \arctan 5 + k\pi\) noktalarıdır. Parantezdeki ifadenin türevi \(-\sin x - \tfrac{1}{5}\cos x\)’tir. \(k\) çiftken \(\sin x\) ve \(\cos x\) pozitif olduğundan bu türev negatiftir. Yani \(f'\) artıdan eksiye geçer ve nokta bir yerel maksimumdur (bkz. Analiz 2). \(k\) tekken aynı akıl yürütme yerel minimum verir. \([0, 20]\) aralığına \(\arctan 5 + 2m\pi\) biçimindeki maksimumlardan \(m = 0, 1, 2\) için olanlar düşer; \(m = 3\) için \(x \approx 20{,}22\) aralığın dışında kalır.
2. adım. Ayrık maksimumları kaydırılmış dilimlerle buluruz. y[1:-1] iç noktaları, y[:-2] sol komşuları, y[2:] sağ komşuları hizalar. (y[1:-1] > y[:-2]) & (y[1:-1] > y[2:]) maskesi iki komşusundan da büyük olan iç noktalarda doğrudur.
3. adım. Maske x[1:-1] noktalarına aittir. Bu yüzden np.flatnonzero ile bulunan indekslere 1 ekleyip x ve y içindeki indeksleri elde ederiz.
import numpy as np
x = np.linspace(0, 20, 2001) # adım 0.01
y = np.exp(-x / 5) * np.sin(x)
peak = (y[1:-1] > y[:-2]) & (y[1:-1] > y[2:]) # iç noktalar için
i = np.flatnonzero(peak) + 1 # x ve y içindeki indeksler
print(x[i])
print(np.round(y[i], 4))
exact = np.arctan(5) + 2 * np.pi * np.arange(3)
print(np.round(exact, 4))
print(np.abs(x[i] - exact).max())Çıktı:
[ 1.37 7.66 13.94]
[0.7451 0.212 0.0604]
[ 1.3734 7.6566 13.9398]
0.003413925875397794Izgara üç yerel maksimum buldu. Bunların kesin noktalardan en büyük uzaklığı yaklaşık \(0{,}0034\)’tür, yani yarım adımdan (\(0{,}005\)) azdır. Maksimum değerler her seferinde aynı oranda küçülür: \(f(x + 2\pi) = e^{-2\pi/5} f(x)\) olduğundan ardışık iki maksimumun oranı \(e^{-2\pi/5} \approx 0{,}285\)’tir. Gerçekten \(0{,}2120 / 0{,}7451 \approx 0{,}285\)’tir. \(\blacksquare\)
Alıştırma 6.7 (Taylor Polinomunu Vandermonde Matrisiyle Hesaplamak) \(e^{-x}\) fonksiyonunun \(0\) etrafındaki üçüncü dereceden Taylor polinomu \(p(x) = \sum_{j=0}^{3} c_j x^j\), \(c_j = (-1)^j/j!\) olsun. \(p\)’yi \(x = 0;\ 0{,}25;\ 0{,}5;\ 0{,}75;\ 1\) noktalarında, \(V_{ij} = x_i^j\) matrisini broadcasting ile kurarak hesaplayın ve \(|p(x) - e^{-x}|\) hatasını \(x^4/24\) sınırıyla karşılaştırın.
Çözüm
x[:, None] (5, 1), k = np.arange(4) ise (4,) şeklindedir. x[:, None]**k ifadesi (5, 4) şekline yayılır ve \(V_{ij} = x_i^j\) kuvvetler tablosunu verir. Bu, NumPy Dizileri bölümünde sütunları yan yana koyarak kurduğumuz Vandermonde matrisidir; broadcasting onu tek ifadeyle verir ve sonucu yine np.vander ile karşılaştırabiliriz. (4,) şekilli \(c\) ile çarpım satırlara yayılır. axis=1 boyunca toplam her nokta için \(\sum_j c_j x_i^j = p(x_i)\) değerini verir.
import math
import numpy as np
x = np.linspace(0, 1, 5) # şekil (5,)
k = np.arange(4) # şekil (4,)
c = np.array([(-1)**j / math.factorial(j) for j in k])
V = x[:, None]**k # şekil (5, 4): V[i, j] = x_i^j
p = (V * c).sum(axis=1) # p(x_i) = sum_j c_j x_i^j
print(V)
print(np.array_equal(V, np.vander(x, 4, increasing=True)))
print(p)
print(np.abs(p - np.exp(-x)))
print(x**4 / 24)Çıktı:
[[1. 0. 0. 0. ]
[1. 0.25 0.0625 0.015625]
[1. 0.5 0.25 0.125 ]
[1. 0.75 0.5625 0.421875]
[1. 1. 1. 1. ]]
True
[1. 0.77864583 0.60416667 0.4609375 0.33333333]
[0. 0.00015495 0.00236399 0.01142905 0.03454611]
[0. 0.00016276 0.00260417 0.01318359 0.04166667]Taylor teoremine göre kalan \(R_3(x) = \frac{f^{(4)}(\xi)}{4!} x^4\)’tür (bkz. Nümerik Analiz). Burada \(f^{(4)}(\xi) = e^{-\xi} \le 1\) olduğundan \(|R_3(x)| \le x^4/24\) olur. Son iki satırı karşılaştırınca her noktada gerçek hatanın bu sınırın altında kaldığı görülür. En büyük hata \(x = 1\)’de, \(0{,}0345 < 0{,}0417\)’dir. \(\blacksquare\)
Alıştırma 6.8 (En Yakın Nokta Çifti) rng = np.random.default_rng(11) üreteciyle birim karede rng.random((200, 2)) noktalarını üretin. Birbirine en yakın iki noktayı uzaklık matrisiyle bulun ve sonucu iç içe iki döngüyle doğrulayın.
Çözüm
Uzaklık matrisini broadcasting ile kurarız. Her çift matriste iki kez geçer (\(D_{ij} = D_{ji}\)) ve köşegen sıfırdır. Bu yüzden köşegeni ve altını np.tril_indices ile np.inf yaparız; geriye yalnız \(i < j\) çiftleri kalır. D.argmin() düzleştirilmiş dizideki indeksi verir, np.unravel_index bunu (satır, sütun) çiftine çevirir.
import numpy as np
rng = np.random.default_rng(11)
P = rng.random((200, 2)) # şekil (200, 2)
diff = P[:, None, :] - P[None, :, :] # şekil (200, 200, 2)
D = np.sqrt((diff**2).sum(axis=-1)) # şekil (200, 200)
D[np.tril_indices(len(P))] = np.inf # i >= j çiftleri dışarıda
i, j = np.unravel_index(D.argmin(), D.shape)
print(i, j, D[i, j])
# Kontrol: iç içe iki döngü, 19900 çift
best, pair = np.inf, None
for a in range(len(P)):
for b in range(a + 1, len(P)):
d = np.hypot(P[a, 0] - P[b, 0], P[a, 1] - P[b, 1])
if d < best:
best, pair = d, (a, b)
print(pair, best == D[i, j])Çıktı:
40 83 0.001200238111042671
(40, 83) TrueEn yakın çift 40 ve 83 numaralı noktalardır; aralarındaki uzaklık yaklaşık \(0{,}0012\)’dir. Döngü aynı çifti ve aynı uzaklığı buldu, ama \(19900\) çifti tek tek dolaşarak. Vektörel çözüm bunun için \(200 \times 200 \times 2\) elemanlı bir ara dizi kurdu. \(\blacksquare\)
Alıştırma 6.9 (Birim Kürenin Hacmi) rng = np.random.default_rng(42) üreteciyle \([-1, 1]^3\) küpünde \(10^6\) düzgün dağılımlı nokta üretin. Birim kürenin içine düşen noktaların oranından kürenin hacmini kestirin ve \(4\pi/3\) ile karşılaştırın.
Çözüm
rng.uniform(-1, 1, size=(n, 3)) her satırı küpte bir nokta olan bir dizi üretir. (P**2).sum(axis=1) <= 1 maskesi noktanın kürenin içinde olup olmadığını söyler. Küpün hacmi \(8\) olduğundan kürenin hacmi yaklaşık \(8 \cdot (\text{içerideki oran})\)’dır.
import numpy as np
rng = np.random.default_rng(42)
n = 10**6
P = rng.uniform(-1, 1, size=(n, 3)) # [-1, 1]^3 küpü, hacmi 8
inside = (P**2).sum(axis=1) <= 1.0
V = 8 * inside.mean()
print(inside.sum(), V)
print(4 * np.pi / 3, abs(V - 4 * np.pi / 3))Çıktı:
523685 4.18948
4.1887902047863905 0.0006897952136091234Bir milyon noktanın \(523685\)’i küreye düştü. Kestirim \(4{,}18948\), gerçek değer \(4\pi/3 \approx 4{,}18879\) ve fark yaklaşık \(7 \cdot 10^{-4}\)’tür. Bir noktanın küreye düşme olasılığı \(p = \pi/6\)’dır. İçerideki noktaların oranı \(\hat{p}\) ise kestirim \(8 \hat{p}\)’dir ve standart sapması \(8\sqrt{p(1 - p)/n} \approx 4 \cdot 10^{-3}\)’tür. Bulunan fark bu beklenen hata düzeyinin rahatça içindedir. \(\blacksquare\)
Alıştırma 6.10 (İkinci Türev İçin Fark Formülü) \([0, 1]\) aralığında \(0{,}01\) adımlı ızgarada \(f(x) = e^x\) değerleri verilsin. İç noktalarda \(f''(x_i) \approx (y_{i+1} - 2y_i + y_{i-1})/h^2\) yaklaşımını dilimlerle hesaplayın. En büyük hatayı \(h^2 e/12\) sınırıyla karşılaştırın.
Çözüm
y[2:], y[1:-1] ve y[:-2] dilimleri sırasıyla \(y_{i+1}\), \(y_i\) ve \(y_{i-1}\) değerlerini iç noktalar için hizalar. Sonuç x[1:-1] noktalarına aittir ve 99 elemanlıdır.
import numpy as np
x = np.linspace(0, 1, 101)
h = x[1] - x[0]
y = np.exp(x)
d2 = (y[2:] - 2 * y[1:-1] + y[:-2]) / h**2 # x[1:-1] noktalarında
err = np.abs(d2 - np.exp(x[1:-1]))
print(len(d2), err.max(), x[1:-1][err.argmax()])
print(h**2 / 12 * np.e)Çıktı:
99 2.242702939581065e-05 0.99
2.2652348570492043e-05Formülün hatası \(\frac{h^2}{12} f^{(4)}(\xi)\)’dir (bkz. Nümerik Analiz). \([0, 1]\) aralığında \(f^{(4)}(x) = e^x \le e\) olduğundan hata en çok \(h^2 e/12 \approx 2{,}27 \cdot 10^{-5}\) olabilir. Bulunan en büyük hata \(2{,}24 \cdot 10^{-5}\)’tir ve \(e^x\)’in en büyük olduğu sağ uca yakın \(x = 0{,}99\) noktasında oluşur. \(\blacksquare\)
Bu bölümde döngüleri dizi işlemlerine çevirmeyi öğrendik. Evrensel fonksiyonlar eleman eleman hesap yapar, axis parametresi indirgemenin yönünü seçer, maskeler ve np.where koşulları döngüsüz işler. Broadcasting ise şekilleri farklı dizileri kopyalamadan birleştirir. Sıradaki NumPy ile Lineer Cebir bölümünde aynı dizileri matris ve vektör olarak ele alacak, matris çarpımı, lineer sistemler ve özdeğerler için NumPy’ın hazır araçlarını kullanacağız.