NumPy: vektorlaşdırma və xətti cəbr

Murad onlayn mağazanın bir milyon satışını emal edir: hər satış üçün məbləği hesablamaq, endirim tətbiq etmək, gündəlik dəyişməni tapmaq lazımdır. for dövrü ilə bu, saniyələr çəkir, NumPy ilə isə millisaniyələr. Sirr vektorlaşdırmadadır: əməliyyatı hər element üçün ayrıca yox, bütöv massiv üçün bir dəfə yazırıq. Bu dərsdə həm də verilənlər elminin riyazi dili olan xətti cəbrlə işləyəcəyik.

Vektorlaşdırma: dövrsüz hesablamalar

Python dövrü niyə yavaşdır? Hər təkrarda interpretator əmri oxuyur, elementin tipini yoxlayır və nəticə üçün yeni float obyekti yaradır. NumPy isə yaddaşda ardıcıl duran eyni tipli ədədlər üzərində bir dəfə tip yoxlaması aparır və prosessorun vektor (SIMD) əmrlərindən istifadə edir. Kodu işə sal və fərqi öz gözünlə gör:

import time
import numpy as np

x = np.random.default_rng(0).random(1_000_000)

t0 = time.perf_counter()
total = 0.0
for v in x:
    total += v * v
t_loop = time.perf_counter() - t0

t0 = time.perf_counter()
total_vec = np.dot(x, x)
t_vec = time.perf_counter() - t0

print(f'loop:       {t_loop:.4f} s')
print(f'vectorized: {t_vec:.4f} s')
print('same result:', bool(np.isclose(total, total_vec)))
print(f'speed-up: about {t_loop / t_vec:.0f} times')

import numpy as np

prices = np.array([12.0, 45.5, 8.9, 120.0, 64.0])
qty = np.array([3, 1, 10, 1, 2])

revenue = prices * qty
discount = np.where(revenue > 100, 0.10, 0.0)
to_pay = (revenue * (1 - discount)).round(2)
print(revenue)
print(to_pay, to_pay.sum())

daily = np.array([520.0, 610.0, 580.0, 700.0, 655.0])
print(np.diff(daily))
print((np.diff(daily) / daily[:-1] * 100).round(1))
print(np.cumsum(daily))

[ 36.   45.5  89.  120.  128. ]
[ 36.   45.5  89.  108.  115.2] 393.7
[ 90. -30. 120. -45.]
[17.3 -4.9 20.7 -6.4]
[ 520. 1130. 1710. 2410. 3065.]

Broadcasting: fərqli formalı massivlər

a * 10 yazanda 10 ədədi sanki a ilə eyni formalı massivə «genişləndirilir». Bu mexanizmə broadcasting (yayılma) deyilir və o, təkcə ədədlərlə deyil, fərqli formalı massivlərlə də işləyir. NumPy formaları sağdan sola müqayisə edir:

1. Ölçülərin sayı azdırsa, formanın soluna 1-lər əlavə olunur: (4,) → (1, 4).
2. Hər ox üzrə iki uzunluq ya bərabər olmalıdır, ya da onlardan biri 1 olmalıdır.
3. Uzunluğu 1 olan ox digər massivin uzunluğuna qədər «uzadılır» (yaddaşda surət çıxarılmadan); nəticənin forması hər ox üzrə maksimumdur.
4. Hər hansı oxda uzunluqlar fərqlidirsə və heç biri 1 deyilsə — ValueError.

Broadcasting-in klassik tətbiqi standartlaşdırmadır: hər sütundan onun ortasını çıxarıb standart meylinə bölürük ki, müxtəlif miqyaslı əlamətlər (boy sm-lə, çəki kq-la) müqayisə oluna bilsin. Forması (4, 2) olan massivdən forması (2,) olan orta çıxılır — qayda 1 və 3 işləyir.

import numpy as np

X = np.array([[170.0, 65.0],
              [160.0, 52.0],
              [182.0, 80.0],
              [175.0, 71.0]])
mean = X.mean(axis=0)
std = X.std(axis=0)
Z = (X - mean) / std
print(mean, std.round(3))
print(Z.round(2))
print(np.allclose(Z.mean(axis=0), 0), Z.std(axis=0).round(6))

[171.75  67.  ] [ 8.012 10.173]
[[-0.22 -0.2 ]
 [-1.47 -1.47]
 [ 1.28  1.28]
 [ 0.41  0.39]]
True [1. 1.]

Sütun vektoru ilə sətir vektorunu birləşdirəndə «hər biri hər biri ilə» cədvəli alınır. a[:, np.newaxis] (və ya a[:, None]) forması (3,) olan massivi (3, 1) sütununa çevirir:

import numpy as np

a = np.arange(1, 4)
b = np.arange(1, 6)
table = a[:, np.newaxis] * b
print(table.shape)
print(table)

M = np.ones((3, 4))
v = np.array([10.0, 20.0, 30.0])
try:
    M + v
except ValueError as e:
    print('Error:', str(e).strip())
print((M + v[:, None])[:, 0])

(3, 5)
[[ 1  2  3  4  5]
 [ 2  4  6  8 10]
 [ 3  6  9 12 15]]
Error: operands could not be broadcast together with shapes (3,4) (3,)
[11. 21. 31.]

Formanı dəyişmək

reshape eyni elementləri başqa formada göstərir; elementlərin sayı dəyişməməlidir, ölçülərdən biri isə -1 ola bilər — NumPy onu özü hesablayır. .T transponə edir, ravel() massivi birölçülü edir (mümkün olanda görünüş qaytarır), np.vstack, np.hstack və np.concatenate isə massivləri yapışdırır.

import numpy as np

a = np.arange(12)
print(a.reshape(3, 4))
print(a.reshape(2, -1).shape)
m = a.reshape(3, 4)
print(m.T.shape, m.ravel()[:5])
print(np.vstack([m[0], m[1]]))
print(np.concatenate([m, m], axis=1).shape)

[[ 0  1  2  3]
 [ 4  5  6  7]
 [ 8  9 10 11]]
(2, 6)
(4, 3) [0 1 2 3 4]
[[0 1 2 3]
 [4 5 6 7]]
(3, 8)

Xətti cəbr: np.linalg

Matrislərin hasili elementbəelement vurma deyil. A * B uyğun elementləri vurur, A @ B (və ya np.dot) isə riyazi matris hasilini hesablayır: A-nın sətri B-nin sütununa skalyar vurulur. Bunun üçün A-nın sütunlarının sayı B-nin sətirlərinin sayına bərabər olmalıdır.

import numpy as np

A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])
print(A * B)
print(A @ B)
print(np.dot(A, B).shape, (np.ones((2, 3)) @ np.ones((3, 5))).shape)

[[ 5 12]
 [21 32]]
[[19 22]
 [43 50]]
(2, 2) (2, 5)

Xətti tənliklər sistemi matris şəklində Ax = b kimi yazılır. det A ≠ 0 olarsa, sistemin yeganə həlli var və A-nın tərs matrisi mövcuddur: x = A⁻¹b. 2 × 2 matris üçün hər şeyi əl ilə hesablamaq olar:

Üç naməlumlu sistemdə (dəftər, qələm, xətkeş) əl ilə hesablama uzanır, NumPy isə onu bir sətirdə həll edir:

import numpy as np

A = np.array([[2.0, 1.0, 3.0],
              [1.0, 3.0, 1.0],
              [4.0, 2.0, 0.0]])
b = np.array([9.5, 5.0, 7.0])

x = np.linalg.solve(A, b)
print(x.round(2))
print(np.allclose(A @ x, b))
print(round(np.linalg.det(A), 2))
x_inv = np.linalg.inv(A) @ b
print(np.allclose(x, x_inv))

[1.5 0.5 2. ]
True
-30.0
True

Məxsusi vektor elə sıfırdan fərqli v vektorudur ki, A matrisi onu yalnız uzadır və ya qısaldır, istiqamətini dəyişmir; uzanma əmsalı λ məxsusi qiymətdir. λ-ları xarakteristik tənlikdən tapırıq. Tətbiqləri: əsas komponentlər üsulu (PCA), Google-un PageRank alqoritmi, körpü və binaların məxsusi rəqs tezlikləri, dinamik sistemlərin dayanıqlığı.

import numpy as np

A = np.array([[2.0, 1.0], [1.0, 2.0]])
values, vectors = np.linalg.eigh(A)
print(values.round(4))
v = vectors[:, 1]
print(np.allclose(A @ v, values[1] * v))
print(round(np.linalg.det(A), 4), round(values.prod(), 4))
print(round(np.trace(A), 4), round(values.sum(), 4))

[1. 3.]
True
3.0 3.0
4.0 4.0

Xətti cəbr maşın öyrənməsinə birbaşa aparır. Nöqtələrdən ən yaxşı düz xətti keçirmək (ən kiçik kvadratlar üsulu) normal tənliyə gətirir. X matrisinə birlərdən ibarət sütun əlavə etsək, sərbəst hədd də avtomatik tapılır:

import numpy as np

rng = np.random.default_rng(7)
x = rng.uniform(0, 10, 50)
y = 2.5 * x + 1.0 + rng.normal(0, 1, 50)
X = np.column_stack([np.ones_like(x), x])
w = np.linalg.solve(X.T @ X, X.T @ y)
print(w.round(3))
w2, *_ = np.linalg.lstsq(X, y, rcond=None)
print(w2.round(3))

[1.236 2.458]
[1.236 2.458]

import numpy as np

X = np.array([[2.0, 50.0, 1.0],
              [4.0, 80.0, 3.0],
              [6.0, 65.0, 2.0],
              [8.0, 95.0, 5.0]])
# min-max scale every column

[[0.   0.   0.  ]
 [0.33 0.67 0.5 ]
 [0.67 0.33 0.25]
 [1.   1.   1.  ]]

import numpy as np

# A: kilograms in each purchase, b: amounts paid (manat)
A = np.array([[2.0, 1.0, 1.0],
              [1.0, 2.0, 1.0],
              [1.0, 1.0, 3.0]])
b = np.array([9.5, 10.0, 13.5])
# solve A x = b and print the prices, then the determinant

[2.  2.5 3. ]
7.0