NumPy：ユニバーサル関数と縮約関数

NumPyのユニバーサル関数（ufunc）¶

NumPyには配列の各要素に同じ処理を高速に適用できる関数群があり，これをユニバーサル関数（ufunc）と呼ぶ．ここでは代表的な ufunc として np.log, np.exp, np.sin, np.cos, np.tan, np.floor, np.isnan を紹介する．

import numpy as np

np.__version__

'2.4.4'

`np.log`¶

np.log は自然対数（底 e の対数）を要素ごとに計算する関数である．入力は正の値を使うのが基本である．

x = np.array([1, np.e, np.e**2, 10])
print("x:", x)
print("np.log(x):", np.log(x))

x: [ 1.          2.71828183  7.3890561  10.        ]
np.log(x): [0.         1.         2.         2.30258509]

`np.exp`¶

np.exp は指数関数 e^x を要素ごとに計算する関数である． np.log と逆の関係を持つため，組み合わせて使われることが多い．

y = np.array([-2, -1, 0, 1, 2])
print("y:", y)
print("np.exp(y):", np.exp(y))
print("np.log(np.exp(y)):", np.log(np.exp(y)))

y: [-2 -1  0  1  2]
np.exp(y): [0.13533528 0.36787944 1.         2.71828183 7.3890561 ]
np.log(np.exp(y)): [-2. -1.  0.  1.  2.]

`np.sin`, `np.cos`, `np.tan`¶

三角関数も ufunc として提供されており，ラジアン単位の角度に対して要素ごとに計算される．角度（度）を扱いたい場合は np.deg2rad でラジアンへ変換してから使うのが分かりやすい．

angles_deg = np.array([0, 30, 45, 60, 90])
angles_rad = np.deg2rad(angles_deg)

print("angles_deg:", angles_deg)
print("angles_rad:", angles_rad)
print("np.sin(angles_rad):", np.sin(angles_rad))
print("np.cos(angles_rad):", np.cos(angles_rad))
print("np.tan(angles_rad):", np.tan(angles_rad))

angles_deg: [ 0 30 45 60 90]
angles_rad: [0.         0.52359878 0.78539816 1.04719755 1.57079633]
np.sin(angles_rad): [0.         0.5        0.70710678 0.8660254  1.        ]
np.cos(angles_rad): [1.00000000e+00 8.66025404e-01 7.07106781e-01 5.00000000e-01
 6.12323400e-17]
np.tan(angles_rad): [0.00000000e+00 5.77350269e-01 1.00000000e+00 1.73205081e+00
 1.63312394e+16]

`np.floor`¶

np.floor は小数を含む値を，値以下の最大の整数へ切り下げる関数である．なお，np.floow という関数は存在せず，正しくは np.floor である．

z = np.array([-1.8, -0.2, 0.0, 1.2, 3.9])
print("z:", z)
print("np.floor(z):", np.floor(z))

z: [-1.8 -0.2  0.   1.2  3.9]
np.floor(z): [-2. -1.  0.  1.  3.]

`np.isnan`¶

np.isnan は要素が NaN（Not a Number）かどうかを判定し，真偽値の配列を返す関数である．欠損値を含むデータの前処理でよく使われる．

w = np.array([0.0, np.nan, 1.5, np.nan, -3.0])
print("w:", w)
print("np.isnan(w):", np.isnan(w))
print("NaNを除いた要素:", w[~np.isnan(w)])

w: [ 0.   nan  1.5  nan -3. ]
np.isnan(w): [False  True False  True False]
NaNを除いた要素: [ 0.   1.5 -3. ]

【問題1】ニューラルネットワークでよく使われるsigmoid関数をNumPyで実装してみよう．¶

\operatorname{sigmoid}(x) = \frac{1}{1+\operatorname{exp}(-x)}

(1)

def sigmoid(x):
    ...

【問題2】ニューラルネットワークでよく使われるsoftmax関数をNumPyで実装してみよう．¶

\mathrm{softmax}(x_i) = \frac{\exp(x_i)}{\sum_{j=1}^{K} \exp(x_j)} \quad (i = 1, \dots, K)

(2)

入力xはベクトルだとしてください．

def softmax(x):
    ...

Reduction（縮約）関数¶

reduction／集約関数：配列全体や軸方向をまとめる関数

`np.sum`¶

X = np.array([
    [2.0, 1.0, 0.1],
    [0.5, 2.5, 1.0]
])

np.sum(X)

np.float64(7.1)

axisについて¶

列ごとに計算，行ごとに計算ができるのがaxis option．

X の形: (2, 3)

[[2.0, 1.0, 0.1],   ← 行 0
 [0.5, 2.5, 1.0]]   ← 行 1
   ↑    ↑    ↑
 列0  列1  列2

axis=0 縦方向（列ごと）に計算

np.sum(X, axis=0) # Xの列の数は3だった．列ごとに合計を計算できるはず．

array([2.5, 3.5, 1.1])

np.sum(X, axis=1) # Xの行の数は2だった．行ごとに合計を計算できるはず．

array([3.1, 4. ])

最小値と最大値（`np.min`, `np.max`）¶

np.min は最小値，np.max は最大値を返す．配列全体だけでなく axis を指定して行ごと・列ごとにも計算できる．

print("np.min(X):", np.min(X))
print("np.max(X):", np.max(X))
print("np.min(X, axis=0):", np.min(X, axis=0))
print("np.max(X, axis=1):", np.max(X, axis=1))

np.min(X): 0.1
np.max(X): 2.5
np.min(X, axis=0): [0.5 1.  0.1]
np.max(X, axis=1): [2.  2.5]

最大値・最小値のインデックス（`np.argmax`, `np.argmin`）¶

np.argmax は最大値をとる位置（インデックス）を返し，np.argmin は最小値をとる位置を返す． axis を指定すると，行ごと・列ごとにインデックスを取得できる．

print("np.argmax(X):", np.argmax(X))
print("np.argmin(X):", np.argmin(X))
print("np.argmax(X, axis=0):", np.argmax(X, axis=0))
print("np.argmin(X, axis=1):", np.argmin(X, axis=1))

np.argmax(X): 5
np.argmin(X): 0
np.argmax(X, axis=0): [1 1 1]
np.argmin(X, axis=1): [0 0]

平均・分散・標準偏差（`np.mean`, `np.var`, `np.std`）¶

np.mean は平均，np.var は分散，np.std は標準偏差を返す．データの中心とばらつきを把握するときに使う基本的な集約関数である．

print("np.mean(X):", np.mean(X))
print("np.var(X):", np.var(X))
print("np.std(X):", np.std(X))
print("np.mean(X, axis=0):", np.mean(X, axis=0))
print("np.std(X, axis=1):", np.std(X, axis=1))

np.mean(X): 1.1833333333333333
np.var(X): 0.6847222222222222
np.std(X): 0.8274794391537607
np.mean(X, axis=0): [1.25 1.75 0.55]
np.std(X, axis=1): [0.77602978 0.84983659]

keepdims=True を付けると？¶

np.sum(X, axis=1, keepdims=True)

array([[3.1],
       [4. ]])

ブロードキャスト¶

ブロードキャストとは形（shape）の異なる配列同士の演算を行うときに、 NumPy が自動的に配列の形をそろえて計算してくれる仕組み

x = np.array([1, 2, 3])
x + 10

array([11, 12, 13])

上の計算では内部で，10 → [10, 10, 10] としてから要素ごとに足している。

では1次元x二次元の配列でも考えてみよう．

X = np.array([
    [1, 2, 3],
    [4, 5, 6]
])

v = np.array([10, 20, 30])

X + v

array([[11, 22, 33],
       [14, 25, 36]])

Xとvは形状が違うが，ちゃんと計算できている．これはブロードキャストのおかげ．

v → [[10, 20, 30],
     [10, 20, 30]]

とみなされてから

X + v → [[11, 22, 33],
          [14, 25, 36]]

こういう計算が行われている．

計算したい二つの配列上がったときに，以下のどちらの条件に当てはまればブロードキャストしてくれる．

次元サイズが同じ
どちらかが 1

「各行に対する値」が列として残るようにできる．配列の形状はプログラムで計算するときに大切なので，これも覚えておこう．

【問題3】多次元配列に対応したsoftmax関数にオーバーフロー対策をつけてみよう．¶

■ softmax関数に入力されるxについて： $$ \mathbf{X}¶

\begin{pmatrix} x_{11} & x_{12} & \cdots & x_{1K} \\ x_{21} & x_{22} & \cdots & x_{2K} \\ \vdots & \vdots & & \vdots \\ x_{N1} & x_{N2} & \cdots & x_{NK} \end{pmatrix}

(3)

\in \mathbb{R}^{N \times K} $$

行：1つのサンプル
列：特徴の数

■ つまりsoftmaxの要素ごとの定義は：

\mathrm{softmax}(\mathbf{X})_{i k} = \frac{\exp(x_{ik})} {\sum_{j=1}^{K} \exp(x_{ij})} \quad (i = 1,\dots,N,\; k = 1,\dots,K)

(4)

各行ごとに softmax が適用される
行ごとに確率の和は 1

■ 数値安定化込みの場合は：

各サンプル $i$ ごとに最大値を引く：

m_i = \max_{1 \le j \le K} x_{ij}

(5)

\mathrm{softmax}(\mathbf{X})_{i k} = \frac{\exp(x_{ik} - m_i)} {\sum_{j=1}^{K} \exp(x_{ij} - m_i)}

(6)

def softmax(x):
    # 入力xはデータ数 $\times$ 特徴数 の配列や，データ数  $\times$ ... $\times$ 特徴数の配列です．xと同じ形状の配列を返してください．
    ...

【問題4】損失関数「二乗和誤差」を実装¶

予測値と正解値の差を二乗して合計し， $\frac{1}{2}$ を掛けた損失を実装してみましょう．回帰でも分類でも，出力と目標のずれを見るときによく使う基本的な指標です．

入力：予測値 $\mathbf{y}$ ，正解値 $\mathbf{t}$ （同じ形状）
出力：スカラーの損失値

E(\mathbf{y},\mathbf{t}) = \frac{1}{2}\sum_{k=1}^{K}(y_k - t_k)^2

(7)

ミニバッチ $N$ 個に対しては，次のように平均しても大丈夫です．

E_{\mathrm{batch}} = \frac{1}{N}\sum_{i=1}^{N}\frac{1}{2}\sum_{k=1}^{K}(y_{ik}-t_{ik})^2

(8)

NumPyでは，配列演算と np.sum を使うとすっきり実装できます．

【問題5】損失関数「クロスエントロピー」を実装¶

分類問題でよく使うクロスエントロピー損失を実装してみましょう．予測確率 $\mathbf{y}$ と正解ラベル（または one-hot） $\mathbf{t}$ から，正解に対する確率が低いほど損失が大きくなります．

入力：予測確率 $\mathbf{y}$ （各要素は 0 以上 1 以下）
入力：正解 $\mathbf{t}$ （one-hot ベクトルを想定）
出力：スカラーの損失値

L(\mathbf{y},\mathbf{t}) = -\sum_{k=1}^{K} t_k\log(y_k)

(9)

数値が不安定にならないように， $\log(0)$ を避ける工夫も入れてみましょう．たとえば $\varepsilon$ を小さい定数として，

L(\mathbf{y},\mathbf{t}) = -\sum_{k=1}^{K} t_k\log(y_k + \varepsilon)

(10)

のように書けます．ミニバッチでは平均損失を返す形にしてみてください．

【問題6】 KL Divergenceを実装¶

2つの確率分布 $\mathbf{p}$ と $\mathbf{q}$ の違いを測る KL Divergence（相対エントロピー）を実装してみましょう． $\mathbf{p}$ を真の分布， $\mathbf{q}$ を近似分布とみなすと，KL Divergence は「 $\mathbf{q}$ が $\mathbf{p}$ をどれくらい近似できているか」を表す量です．

入力：確率分布 $\mathbf{p}, \mathbf{q}$ （同じ形状，和は1）
出力：スカラーの非負値（理論上 $\ge 0$ ）

D_{\mathrm{KL}}(\mathbf{p}\|\mathbf{q}) = \sum_{k=1}^{K} p_k\log\frac{p_k}{q_k}

(11)

ミニバッチ $N$ 個の場合は，次のように平均してもよいです．

D_{\mathrm{KL}}^{\mathrm{batch}} = \frac{1}{N}\sum_{i=1}^{N}\sum_{k=1}^{K} p_{ik}\log\frac{p_{ik}}{q_{ik}}

(12)

ゼロ割りや $\log(0)$ を避けるために，小さい $\varepsilon$ を加える実装にしてみてください．

【問題7】Sigmoidクラスを実装¶

データ数 $\times$ 特徴数の配列 X を受け取って，sigmoid の計算を行うメソッド forward と，合成関数の微分を用いて入力 X 側への変化率を計算するメソッド backward を持ったクラス Sigmoid を作成してください．

forward(X) は，各要素に sigmoid 関数を適用した配列を返す．
backward(dout) は，出力側から与えられる配列 dout を受け取り，合成関数の微分を使って入力 X に関する変化率を計算して返す．

合成関数の流れは次のように表せる：

A \xrightarrow{\text{別の関数}} X \xrightarrow{\mathrm{sigmoid}} Y \xrightarrow{\text{別の関数}} D

(13)

ここで，Sigmoid クラスは $X \xrightarrow{\mathrm{sigmoid}} Y$ の部分を担当する．

sigmoid 関数は次式で与えられる．

\operatorname{sigmoid}(x) = \frac{1}{1 + \exp(-x)}

(14)

また， $y = \operatorname{sigmoid}(x)$ とすると，

\frac{dy}{dx} = y(1-y)

(15)

が成り立つ．（本当にこうなるのか，紙と鉛筆を出して導出してみよう）

したがって，合成関数の微分より backward(dout) では各要素について

dX = dout \cdot y(1-y)

(16)

を計算すればよい．必要であれば，forward で計算した出力をクラス内に保存して backward で利用してよい．

class Sigmoid:
    def __init__(self):
        self._output = None

    def forward(self, x):
        ...
        return self._output
    
    def backward(self, dout):
        ...
        return dx

# テストコード

x = np.array([[0.0]])
dout = np.array([[1.0]])

layer = Sigmoid()

y = layer.forward(x)
dx = layer.backward(dout)

print("y =", y)
print("dx =", dx)

assert np.allclose(y, np.array([[0.5]]))
assert np.allclose(dx, np.array([[0.25]]))

print("1x1 test passed")

# 追加テストコード（2×2）
x = np.array([
    [0.0, 1.0],
    [-1.0, 2.0],
])
dout = np.array([
    [1.0, 2.0],
    [3.0, 4.0],
])

layer = Sigmoid()

y = layer.forward(x)
dx = layer.backward(dout)

expected_y = np.array([
    [0.5, 0.73105858],
    [0.26894142, 0.88079708],
])
expected_dx = np.array([
    [0.25, 0.39322387],
    [0.5898358, 0.41997434],
])

print("y =")
print(y)
print("dx =")
print(dx)

assert np.allclose(y, expected_y)
assert np.allclose(dx, expected_dx)

print("2x2 test passed")

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[8], line 9
      5 
      6 layer = Sigmoid()
      7 
      8 y = layer.forward(x)
----> 9 dx = layer.backward(dout)
     10 
     11 print("y =", y)
     12 print("dx =", dx)

Cell In[7], line 11, in Sigmoid.backward(self, dout)
      9     def backward(self, dout):
     10         ...
---> 11         return dx

NameError: name 'dx' is not defined

【発展問題1】線形回帰を実装¶

特徴量 $\mathbf{x}$ から連続値 $\hat{y}$ を予測する線形回帰を NumPy で実装してみましょう．まずは予測式を作って，そのあと二乗和誤差が小さくなるようにパラメータを更新していきます．

モデル¶

1サンプル $\mathbf{x}\in\mathbb{R}^{D}$ に対して，

\hat{y} = \mathbf{w}^{\top}\mathbf{x} + b

(17)

ミニバッチ $\mathbf{X}\in\mathbb{R}^{N\times D}$ では，

\hat{\mathbf{y}} = \mathbf{X}\mathbf{w} + b

(18)

損失（二乗和誤差の平均）¶

J(\mathbf{w}, b) = \frac{1}{2N}\sum_{i=1}^{N}(\hat{y}_i - y_i)^2

(19)

勾配¶

\frac{\partial J}{\partial \mathbf{w}} = \frac{1}{N}\mathbf{X}^{\top}(\hat{\mathbf{y}}-\mathbf{y}), \qquad \frac{\partial J}{\partial b} = \frac{1}{N}\sum_{i=1}^{N}(\hat{y}_i-y_i)

(20)

学習率 $\eta$ を使うと，更新は次のように書けます．

\mathbf{w}\leftarrow\mathbf{w}-\eta\frac{\partial J}{\partial \mathbf{w}}, \qquad b\leftarrow b-\eta\frac{\partial J}{\partial b}

(21)

この更新を何度も行うことでより良いパラメータを見つけます．

この式をもとに，LinearRegressionクラスを作ってください．ただし，以下の条件を守ってください．

式中の $x$ , $y$ に相当する配列をfitメソッドの引数X,yとして受け取り，このメソッド中で学習によって獲得したパラメータをself._W, self._bに保存します．
predictメソッドで数式中の $x$ に相当する配列を引数Xとして受け取り，予測値を返してください．
パラメータの学習には以下のコードで生成するX_train,y_trainを利用して，predictではX_test, y_testを利用してください．
予測がy_testとどれだけ違うかを二乗和誤差で算出してください。
作成したクラスを使って回帰直線を作成し，fitで利用したデータの散布図上に回帰直線を描いたグラフを作って下さい．

from jaxtyping import Float
from typing import Final

def make_linear_regression_data_2d(
    n_samples: int = 200,
    noise_std: float = 0.3,
    seed: int = 42,
) -> tuple[Float[np.ndarray, "n 2"], Float[np.ndarray, "n"]]:
    """Generate a 2D synthetic dataset for linear regression.

    Args:
        n_samples: Number of samples.
        noise_std: Standard deviation of Gaussian noise.
        seed: Random seed.

    Returns:
        A tuple of feature matrix X and target vector y.
        X has shape (n_samples, 2), y has shape (n_samples,).
    """
    rng = np.random.default_rng(seed)

    X: Float[np.ndarray, "n 2"] = rng.uniform(-2.0, 2.0, size=(n_samples, 2))
    true_w: Final[np.ndarray] = np.array([2.5, -1.2], dtype=np.float64)
    true_b: Final[float] = 0.7
    noise: Float[np.ndarray, "n"] = rng.normal(0.0, noise_std, size=n_samples)

    y: Float[np.ndarray, "n"] = X @ true_w + true_b + noise
    return X, y

def train_test_split_regression_data(
    X: Float[np.ndarray, "n d"],
    y: Float[np.ndarray, "n"],
    test_size: float = 0.2,
    seed: int = 42,
) -> tuple[
    Float[np.ndarray, "n_train d"],
    Float[np.ndarray, "n_test d"],
    Float[np.ndarray, "n_train"],
    Float[np.ndarray, "n_test"],
]:
    """Split regression data into train and test sets.

    Args:
        X: Feature matrix of shape (n_samples, n_features).
        y: Target vector of shape (n_samples,).
        test_size: Fraction of samples used for the test set.
        seed: Random seed.

    Returns:
        A tuple of X_train, X_test, y_train, y_test.
    """
    rng = np.random.default_rng(seed)
    n_samples = X.shape[0]
    indices = rng.permutation(n_samples)

    n_test = int(n_samples * test_size)
    test_indices = indices[:n_test]
    train_indices = indices[n_test:]

    X_train = X[train_indices]
    X_test = X[test_indices]
    y_train = y[train_indices]
    y_test = y[test_indices]

    return X_train, X_test, y_train, y_test


X_all, y_all = make_linear_regression_data_2d()

X_train, X_test, y_train, y_test = train_test_split_regression_data(
    X_all,
    y_all,
    test_size=0.2,
    seed=42,
)

print("X_train:", X_train.shape)
print("X_test:", X_test.shape)
print("y_train:", y_train.shape)
print("y_test:", y_test.shape)

X_train: (160, 2)
X_test: (40, 2)
y_train: (160,)
y_test: (40,)

class LinearRegression:
    def __init__(self, input_dim: int, learning_rate: float = 0.01, num_iterations: int = 1000):
        self._W = np.random.rand(input_dim) # 初期化
        self._b = float(np.random.rand(1)) # 初期化
        self.learning_rate = learning_rate
        self.num_iterations = num_iterations

    def fit(self, X, y):
        # Xはデータ数 $\times$ 特徴数の配列で，yはデータ数の配列です．重みとバイアスを更新してください．
        ...
        return self
    
    def predict(self, X):
        # Xはデータ数 $\times$ 特徴数の配列です．予測値を返してください．
        return ...

作成したクラスを使ってみよう

散布図と回帰直線が描かれたグラフを作成してみよう．matplotlibの使い方は自分で調べること．

NumPyのユニバーサル関数（ufunc）¶

np.log¶

np.exp¶

np.sin, np.cos, np.tan¶

np.floor¶

np.isnan¶