Pythonで学ぶ数値計算のアルゴリズムと実装第8回：スパース性を活かす：帯行列の高速計算アルゴリズム

こんにちは、JS2IIUです。
数値計算の連載も後半戦に入りました。前回までの記事では、密行列（要素のほとんどが0ではない行列）を対象としたLU分解やコレスキー分解を扱ってきました。これらは非常に汎用性が高い手法ですが、計算量が $O(n^3)$ であり、行列のサイズ $n$ が大きくなると計算コストが爆発的に増大するという課題がありました。

しかし、実世界の多くの工学問題や物理シミュレーション、あるいは特定の構造を持つニューラルネットワーク（畳み込み演算の数学的背景など）において、現れる行列のほとんどの要素が「0」であるケースが多々あります。このような行列を「スパース行列（疎行列）」と呼びます。

今回は、スパース行列の中でも特に重要な「帯行列（Band Matrix）」に焦点を当てます。帯行列とは、非ゼロ要素が主対角線の周囲にのみ集中している行列のことです。この構造を賢く利用することで、計算量とメモリ使用量を劇的に削減する方法を、PythonとPyTorchを通じて学んでいきましょう。今回もよろしくお願いします。

1. はじめに：なぜ帯行列が重要なのか？
2. 帯行列の定義と構造
3. メモリ節約の鍵：帯行列の格納方式
4. 帯LU分解アルゴリズムの原理
1. 帯LU分解の計算手順
  1. 数式イメージ
5. PyTorchによる実践：帯行列計算のシミュレーション
1. コードの解説
6. 帯コレスキー分解：対称帯行列のさらなる高速化
7. 帯行列の逆行列に関する注意点
1. 教訓：逆行列を求めてはいけない
8. 大規模計算におけるメモリ管理の重要性
まとめ
連載記事、各回へのリンク

1. はじめに：なぜ帯行列が重要なのか？

大規模な数値計算やシミュレーションの現場では、行列のサイズが1万次元、10万次元と非常に大きくなることが珍しくありません。第6回や第7回で学んだ通り、こうした「密行列」（ほとんどの要素が非ゼロ）のLU分解を行う場合、計算量は $O(n^3)$ にもなり、例えば10,000次元の行列では約3,300億回もの演算が必要となります。これは現代のコンピュータでも膨大な負荷です。

しかし、現実の多くの工学・物理・AIの問題では、行列のほとんどの要素が「0」であることが多いのです。特に、主対角線の近くにだけ非ゼロ要素が集中し、それ以外はすべて0となる「帯行列（バンド行列）」が頻繁に現れます。

帯行列が重要となる代表的な例は次の通りです。

微分方程式の数値解法: たとえば熱伝導や波動方程式を差分法で離散化すると、各格子点は隣接点としか相互作用しないため、行列は主対角線の周囲だけに非ゼロ要素を持つ帯状の構造になります。
スプライン補間: データの滑らかな補間を行う際に現れる方程式も、隣接点のみが関係するため帯行列となります。
信号処理: 1次元信号のフィルタリングや畳み込み演算も、行列で表現すると帯行列の形になります。

このような場合、行列の大部分を占める「0」を無視し、非ゼロの「帯」部分だけを格納・計算することで、計算量を劇的に削減できます。たとえば帯幅$w$の帯行列なら、計算量は $O(n^3)$ から $O(n \cdot w^2)$ へと大幅に減少します（$w$は帯の幅）。

つまり、帯行列の構造を活かすことは「計算を速くする」だけでなく、「そもそも大規模な問題を現実的なメモリ・時間で解くための必須テクニック」なのです。スパース性（疎性）を見抜き、帯行列として扱うことが、現代の数値計算・AI実装のスケーラビリティを大きく左右します。

2. 帯行列の定義と構造

帯行列（バンド行列）とは、主対角線の近くにだけ非ゼロ要素が集中し、それ以外の要素がすべて0であるような特殊な構造を持つ行列です。

より厳密には、$n \times n$ 行列 $A$ について、主対角線から下に $k_1$ 本、上に $k_2$ 本離れた範囲外の要素がすべて $0$ であるとき、$A$ は帯行列と呼ばれます。

$$A_{ij} = 0 \quad (\text{if } j < i – k_1 \text{ or } j > i + k_2)$$

ここで：

$k_1$：下帯幅（Lower bandwidth）…主対角線から下に何本分まで非ゼロ要素が許されるか
$k_2$：上帯幅（Upper bandwidth）…主対角線から上に何本分まで非ゼロ要素が許されるか
$w = k_1 + k_2 + 1$：全帯幅（Total bandwidth）…非ゼロが存在しうる横幅

この構造を図でイメージすると、行列の中央（主対角線）を中心に、上下$k_1$本・$k_2$本だけ太い帯が走っており、それ以外はすべて0です。

例えば、$k_1=1, k_2=1$ の場合は「三重対角行列」と呼ばれ、各行に最大3つしか非ゼロ要素がありません：

$$
\begin{bmatrix}
* & * & 0 & 0 & \cdots & 0 \\
* & * & * & 0 & \cdots & 0 \\
0 & * & * & * & \cdots & 0 \\
\vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\
0 & 0 & 0 & * & * & * \\
0 & 0 & 0 & 0 & * & *
\end{bmatrix}
$$

このように、帯幅が小さいほど、計算やメモリ消費が大幅に削減できます。たとえば$1000 \times 1000$の三重対角行列なら、全体で100万要素のうち、実際に計算に関わるのは約3,000個だけです。

3. メモリ節約の鍵：帯行列の格納方式

帯行列の最大の利点のひとつは、メモリ消費を劇的に削減できる点です。通常の $N \times N$ 行列は $N^2$ 個の要素をすべて格納しますが、帯行列の場合、非ゼロ要素は主対角線の周囲の「帯」部分だけに集中しています。

このため、帯行列は「帯格納形式（Banded Storage）」と呼ばれる特殊な方法で保存するのが一般的です。具体的には、$A$ の全要素を保存するのではなく、非ゼロとなりうる帯部分だけを $(k_1 + k_2 + 1) \times N$ の矩形配列として格納します。

例えば、$k_1=1, k_2=1$ の三重対角行列（tridiagonal matrix）なら、

$$
A = \begin{bmatrix}
a_{11} & a_{12} & 0 & \cdots & 0 \\
a_{21} & a_{22} & a_{23} & \cdots & 0 \\
0 & a_{32} & a_{33} & \cdots & 0 \\
\vdots & \vdots & \vdots & \ddots & \vdots \\
0 & 0 & 0 & \cdots & a_{nn}
\end{bmatrix}
$$

このとき帯格納形式では、

$$
B = \begin{bmatrix}
ext{上帯} & a_{12} & a_{23} & \cdots & a_{n-1,n} \\
ext{主対角} & a_{11} & a_{22} & \cdots & a_{nn} \\
ext{下帯} & a_{21} & a_{32} & \cdots & a_{n, n-1}
\end{bmatrix}
$$

のように、各列ごとに「上帯」「主対角」「下帯」の成分だけを $3 \times N$ の配列として保存します。

一般の帯行列（帯幅 $w = k_1 + k_2 + 1$）では、

$$
B_{ij} = A_{i + j – k_1, j} \quad (0 \leq i < w,\ 0 \leq j < N)
$$

といった対応で、$A$ の帯部分だけを $w \times N$ の配列 $B$ に詰めていきます。

この方法により、$N$ が大きくても帯幅 $w$ が小さければ、必要なメモリ量は $O(Nw)$ で済みます。たとえば $N=10^5$, $w=5$ なら、$10^5 \times 5 = 50$ 万要素だけで済み、$N^2=10^{10}$ 要素を保存する密行列に比べて圧倒的な節約となります。

Pythonの SciPy ライブラリ（scipy.linalg.solve_banded など）では、この帯格納形式が標準的に使われています。PyTorchには現状、帯行列専用のストレージ型はありませんが、自作アルゴリズムでもこの形式を意識することで、GPUメモリの消費を数百分の一に抑えることができます。

4. 帯LU分解アルゴリズムの原理

帯行列に対してLU分解を行うと、非常に大きなメリットがあります。それは、「分解後の $L$（下三角行列）と $U$（上三角行列）も、元の $A$ と同じ帯幅（またはそれに近い幅）を維持する」という性質です。

通常のLU分解では、ガウス消去の過程で本来ゼロだった場所に新たな非ゼロ要素が生まれる「フィルイン（Fill-in）」が問題となり、スパース性が失われてしまうことがあります。しかし帯行列の場合、消去操作の影響範囲が帯の中に限定されるため、$L$と$U$も同じ帯幅を持つ（またはごくわずかしか広がらない）という嬉しい性質があります。

この理由は、帯行列の各行・各列で非ゼロ要素が存在する範囲があらかじめ決まっており、消去操作（行の引き算）もその範囲内でしか影響を及ぼさないからです。

帯LU分解の計算手順

各列 $k$ について、主対角線から下に $k_1$ 行分だけ（$i = k+1$ から $i = \min(k+k_1, n-1)$ まで）を消去対象とします。
各行の更新も、主対角線から右に $k_2$ 列分だけ（$j = k+1$ から $j = \min(k+k_2, n-1)$ まで）しか計算しません。

この制約を加えることで、計算量は$O(n \cdot w^2)$（$w$は帯幅）に抑えられ、$O(n^3)$の密行列LU分解に比べて圧倒的に効率的です。

数式イメージ

例えば三重対角行列（$k_1 = k_2 = 1$）の場合、$k$列目の消去操作は$k+1$行目だけに影響し、$L$と$U$の非ゼロパターンも三重対角のまま保たれます。

一般の帯行列でも、
$$
L_{ik} = \frac{A_{ik}^{(k)}}{A_{kk}^{(k)}} \quad (i = k+1, \ldots, k+k_1)
$$
$$
A_{ij}^{(k+1)} = A_{ij}^{(k)} – L_{ik} A_{kj}^{(k)} \quad (j = k+1, \ldots, k+k_2)
$$
のように、更新範囲が帯の中に限定されます。

このようにして、帯行列のLU分解は「スパース性を保ったまま」効率的に進行し、計算コスト・メモリ消費の両面で大きな恩恵をもたらします。

5. PyTorchによる実践：帯行列計算のシミュレーション

ここでは、帯行列の構造がいかに計算を効率化できるかを理解するために、PyTorchを用いて「帯幅を限定したLU分解」のロジックを実装してみましょう。

以下のコードでは、通常の全行列を対象とした処理と、帯の範囲内のみを計算する処理を対比させています。

Python

import torch
import time

def band_lu_decomposition(A, k1, k2):
    """
    帯行列 A に対して、帯幅を考慮した LU 分解を行う (簡易版)
    A: n x n の正方行列 (帯行列であることを前提とする)
    k1: 下帯幅
    k2: 上帯幅
    """
    n = A.shape[0]
    # 元の行列を破壊しないようにコピー
    LU = A.clone().float()

    # 計算開始
    start_time = time.time()

    for k in range(n):
        # 消去範囲を限定: 下に k1 行分だけ
        lower_limit = min(k + k1 + 1, n)
        for i in range(k + 1, lower_limit):
            # 帯行列の性質上、LU[k, k] が 0 になることは少ない (ピボットなしの場合)
            factor = LU[i, k] / LU[k, k]
            LU[i, k] = factor

            # 更新範囲を限定: 右に k2 列分だけ
            upper_limit = min(k + k2 + 1, n)
            LU[i, k+1:upper_limit] -= factor * LU[k, k+1:upper_limit]

    end_time = time.time()
    return LU, end_time - start_time

# テスト用データの作成 (n=1000, 帯幅3の三項対角行列)
n = 1000
k1, k2 = 1, 1
A_band = torch.zeros((n, n))

# 対角成分と上下1本の帯を埋める
for i in range(n):
    A_band[i, i] = 4.0
    if i > 0: A_band[i, i-1] = 1.0
    if i < n-1: A_band[i, i+1] = 1.0

# 1. 通常の(密行列を想定した)計算 (比較用としてk1, k2を最大値に設定)
_, dense_time = band_lu_decomposition(A_band, n, n)
print(f"Dense LU processing time: {dense_time:.4f} seconds")

# 2. 帯構造を活かした計算
_, band_time = band_lu_decomposition(A_band, k1, k2)
print(f"Band LU processing time:  {band_time:.4f} seconds")

print(f"Speedup: {dense_time / band_time:.2f}x")

import torch
import time

def band_lu_decomposition(A, k1, k2):
    """
    帯行列 A に対して、帯幅を考慮した LU 分解を行う (簡易版)
    A: n x n の正方行列 (帯行列であることを前提とする)
    k1: 下帯幅
    k2: 上帯幅
    """
    n = A.shape[0]
    # 元の行列を破壊しないようにコピー
    LU = A.clone().float()

    # 計算開始
    start_time = time.time()

    for k in range(n):
        # 消去範囲を限定: 下に k1 行分だけ
        lower_limit = min(k + k1 + 1, n)
        for i in range(k + 1, lower_limit):
            # 帯行列の性質上、LU[k, k] が 0 になることは少ない (ピボットなしの場合)
            factor = LU[i, k] / LU[k, k]
            LU[i, k] = factor

            # 更新範囲を限定: 右に k2 列分だけ
            upper_limit = min(k + k2 + 1, n)
            LU[i, k+1:upper_limit] -= factor * LU[k, k+1:upper_limit]

    end_time = time.time()
    return LU, end_time - start_time

# テスト用データの作成 (n=1000, 帯幅3の三項対角行列)
n = 1000
k1, k2 = 1, 1
A_band = torch.zeros((n, n))

# 対角成分と上下1本の帯を埋める
for i in range(n):
    A_band[i, i] = 4.0
    if i > 0: A_band[i, i-1] = 1.0
    if i < n-1: A_band[i, i+1] = 1.0

# 1. 通常の(密行列を想定した)計算 (比較用としてk1, k2を最大値に設定)
_, dense_time = band_lu_decomposition(A_band, n, n)
print(f"Dense LU processing time: {dense_time:.4f} seconds")

# 2. 帯構造を活かした計算
_, band_time = band_lu_decomposition(A_band, k1, k2)
print(f"Band LU processing time:  {band_time:.4f} seconds")

print(f"Speedup: {dense_time / band_time:.2f}x")

Plaintext

Dense LU processing time: 23.8431 seconds
Band LU processing time:  0.0424 seconds
Speedup: 562.74x

Dense LU processing time: 23.8431 seconds
Band LU processing time:  0.0424 seconds
Speedup: 562.74x

コードの解説

範囲限定ループ: lower_limit と upper_limit を計算することで、行列全体の $n$ 回ループではなく、帯幅分（$k_1, k_2$）のループに抑えています。
スライシングの活用: LU[i, k+1:upper_limit] のようにPyTorchのスライシングを使うことで、GPUやCPUのベクトル演算能力を引き出しつつ、不要な0へのアクセスを遮断しています。
計算速度の差: 行列サイズ $n$ が大きくなればなるほど、この速度差は顕著になります。

6. 帯コレスキー分解：対称帯行列のさらなる高速化

行列が「対称（$A = A^T$）」かつ「帯行列」である場合、第7回で学んだコレスキー分解の帯行列版、「帯コレスキー分解（Banded Cholesky Decomposition）」が極めて強力な手法となります。

帯コレスキー分解は、対称正定値な帯行列 $A$ を、
$$A = LL^T$$
という形で、$L$ を同じ帯幅を持つ下三角行列として分解する方法です。

帯コレスキー分解の計算量とメモリ効率

対称帯行列では $k_1 = k_2 = k$ となり、$A$ の非ゼロ要素は主対角線の上下 $k$ 本の帯に集中しています。$L$ も同じ帯幅 $k$ を持つ下三角帯行列となります。

この場合、計算量は $O(n \cdot k^2)$ となり、密行列の $O(n^3)$ に比べて圧倒的に高速です。必要なメモリも $O(nk)$ で済みます。

帯コレスキー分解の計算手順

帯コレスキー分解では、$L$ の成分は次のように計算されます：

対角成分：
$$l_{jj} = \sqrt{a_{jj} – \sum_{s=1}^{k} l_{js}^2}$$
（ただし $l_{js}$ は $a_{js}$ が非ゼロとなる帯内の成分のみ）
非対角成分（$i > j$ かつ $|i-j| \leq k$）：
$$l_{ij} = \frac{1}{l_{jj}} \left(a_{ij} – \sum_{s=1}^{j-1} l_{is}l_{js}\right)$$
（この和も帯内の成分だけを計算すればよい）

このように、計算対象が帯の範囲に限定されるため、ループ回数・演算量ともに大幅に削減されます。

実用例と応用分野

帯コレスキー分解は、

構造力学における大規模トラス構造の有限要素法解析
音響信号処理や自己回帰モデルの自己相関行列
微分方程式の差分法による離散化
など、対称正定値な帯行列が現れる多くの分野で「ゴールドスタンダード」として使われています。

この手法を使うことで、従来は計算・メモリ的に不可能だった大規模問題も、現実的な時間とリソースで解くことが可能になります。

7. 帯行列の逆行列に関する注意点

ここで、数値計算の中級者が陥りやすい重要な罠についてお話しします。
「帯行列の逆行列は、帯行列ではない（通常は密行列になる）」という事実です。

行列 $A$ がどんなに細い帯行列であっても、その逆行列 $A^{-1}$ を計算すると、ほぼすべての要素が非ゼロで埋まってしまいます。

教訓：逆行列を求めてはいけない

大規模な帯行列を扱う際、逆行列を陽に計算しようとすると、スパース性を失いメモリが溢れます。

「逆行列を求めてからベクトル $b$ を掛ける」のではなく、
「帯LU分解（または帯コレスキー分解）を行ってから代入法で解く」

この原則を守ることが、大規模数値計算を成功させる鉄則です。

8. 大規模計算におけるメモリ管理の重要性

現代のAI開発や科学技術計算において、GPUやメインメモリは非常に貴重なリソースです。特にディープラーニングや大規模シミュレーションでは、扱う行列やテンソルのサイズが飛躍的に大きくなっています。

例えば、$n = 10^6$（100万次元）のベクトルや行列を扱う場合、密行列として $n \times n = 10^{12}$ 要素を保存しようとすると、単精度（4バイト）でも約4TBものメモリが必要です。これは最新のGPUや多くのサーバーでも到底実現できません。

しかし、現実の多くの問題では「帯行列」や「スパース行列」として扱える場合が多く、帯幅 $w = 100$ であれば、必要なメモリは $n \times w = 10^8$ 要素、約400MB程度で済みます。これはノートPCや一般的なGPUでも十分扱えるサイズです。

このように、帯行列アルゴリズムを活用することで、

これまで不可能だった大規模な物理シミュレーションやAIモデルの学習
超高解像度データや長大な時系列データの処理
科学技術計算・工学設計の現場でのリアルタイム解析
などが現実的なリソースで実現できるようになります。

また、メモリ消費が減ることで、

より多くのデータやモデルを同時に扱える
バッチサイズやシミュレーション規模を拡大できる
計算速度も向上し、開発サイクルが短縮される
といった副次的なメリットも得られます。

帯行列・スパース行列の考え方は、単なる「高速化テクニック」ではなく、「大規模データ時代の数値計算・AI実装の根幹」となる重要な知識です。メモリ管理の観点からも、ぜひ意識して活用していきましょう。

まとめ

第8回では、行列の「形」を利用した最適化手法である帯行列アルゴリズムについて学びました。

帯行列は、対角線周囲にのみ要素を持つ、構造化されたスパース行列である。
帯格納形式を用いることで、メモリ消費を劇的に抑えられる。
帯LU/コレスキー分解は、計算量を $O(n^3)$ から $O(n \cdot w^2)$ に短縮する。
フィルインの抑制こそが、帯行列アルゴリズムの真骨頂である。
逆行列の密化を避け、常に分解と代入で解を求めるべきである。

スパース性を理解し、活用できるようになると、あなたの書くPythonコードの「扱えるデータ規模」が一気に広がります。

次回は、データの当てはめ（フィッティング）に欠かせない「最小二乗法」と、それを数値的に極めて安定して解くための「ハウスホルダー変換」および「QR分解」について解説します。統計学やAIの根幹を支える、非常に強力な手法です。お楽しみに。

最後まで読んでいただきありがとうございます。

連載記事、各回へのリンク

Pythonで学ぶ数値計算のアルゴリズムと実装第1回：多項式計算とテイラー展開 | アマチュア無線局JS2IIU
https://js2iiu.com/2025/12/31/python_numerical_01/

Pythonで学ぶ数値計算のアルゴリズムと実装第2回：多項式補間とチェビシェフ補間の戦略 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/01/python_numerical_02/

Pythonで学ぶ数値計算のアルゴリズムと実装第3回：積分の基礎：中点則、台形則からシンプソン則まで | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/01/python_numerical_03/

Pythonで学ぶ数値計算のアルゴリズムと実装第4回：高度な数値積分：二重指数型（DE）公式と広義積分 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/01/python_numerical_04/

Pythonで学ぶ数値計算のアルゴリズムと実装第5回：計算の「正しさ」を評価する：ノルムと条件数 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/01/python_numeric_05/

Pythonで学ぶ数値計算のアルゴリズムと実装第6回：連立一次方程式の直截解法：ガウス消去法とピボット選択 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/01/python_numerical_06/

Pythonで学ぶ数値計算のアルゴリズムと実装第8回：スパース性を活かす：帯行列の高速計算アルゴリズム | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/02/python_numerical_08/

Pythonで学ぶ数値計算のアルゴリズムと実装第9回：最小二乗法とQR分解：ハウスホルダー変換による安定化 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/02/python_numerial_09/

Pythonで学ぶ数値計算のアルゴリズムと実装第10回：非線形方程式の探索：2分法からニュートン法まで | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/03/python_numerical_10/

Pythonで学ぶ数値計算のアルゴリズムと実装第11回：行列の本質を捉える：固有値問題の数値解法 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/03/python_numerical_11/

Pythonで学ぶ数値計算のアルゴリズムと実装第12回：現象の変化を追う：ルンゲ・クッタ法による常微分方程式の解法 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/03/python_numerical_12/

Pythonで学ぶ数値計算のアルゴリズムと実装第13回：空間の熱と振動を解く：偏微分方程式の差分近似と反復解法 | アマチュア無線局JS2IIU
https://js2iiu.com/2026/01/03/python_numerical_13/

月	火	水	木	金	土	日
			1	2	3	4
5	6	7	8	9	10	11
12	13	14	15	16	17	18
19	20	21	22	23	24	25
26	27	28	29	30	31