第10章: 区間演算

Interval Arithmetic

歴史と動機

数値計算において、浮動小数点演算は各ステップで丸め誤差を生じる。 例えば IEEE 754 倍精度では約 $10^{-16}$ の相対誤差が毎回発生し、 数千ステップの計算を経ると誤差がどこまで蓄積するかは自明でない。 従来はこの問題に対し、手作業による誤差解析が行われてきた。 すなわち、各演算の丸め誤差を個別に評価し、それらを累積して最終誤差の上界を導出するという手法である。 しかし、この方法は煩雑で人為的ミスが生じやすく、 大規模な数値計算に対しては非現実的であった。

Ramon Moore と区間演算の誕生

1960年代、Ramon E. Moore は博士論文および著書 "Interval Analysis"（1966年）において、 区間演算（interval arithmetic）を体系的に提唱した。 核となるアイデアは単純である： 単一の浮動小数点数の代わりに、真の値を含むことが保証された区間を用いて計算を行う。 演算のたびに区間が正しく伝搬されるよう定義すれば、 最終結果の区間が自動的に誤差の上界を与える。

区間演算の動機をまとめると以下のようになる：

  • 自動化：手作業の誤差解析を不要にし、計算機自身が誤差範囲を追跡する
  • 厳密性：得られた区間は真の値を必ず含むことが数学的に保証される
  • 汎用性：四則演算だけでなく、初等関数や方程式の求解にも拡張可能である
  • 検証：数値計算の結果に対し、数学的に厳密な存在証明や不在証明を与えることができる

Moore 以降、区間演算は精力的に発展し、 Alefeld & Herzberger (1983)、Neumaier (1990)、Tucker (2011) など多くの教科書が著された。 現在では、計算機援用証明（computer-assisted proof）の基盤技術として、 微分方程式の解の存在証明やカオス力学系の解析など幅広い分野で活用されている。

区間演算の定義

基本概念

定義（閉区間）

実数 $\underline{x}, \overline{x}$ に対し $\underline{x} \leq \overline{x}$ のとき、

$$\mathbf{x} = [\underline{x}, \overline{x}] = \{ x \in \mathbb{R} \mid \underline{x} \leq x \leq \overline{x} \}$$

を閉区間と呼ぶ。$\underline{x}$ を下端（infimum）、$\overline{x}$ を上端（supremum）という。 すべての閉区間の集合を $\mathbb{IR}$ と書く。

定義（区間の諸量）

区間 $\mathbf{x} = [\underline{x}, \overline{x}]$ に対し、以下を定義する：

  • 幅（width）：$w(\mathbf{x}) = \overline{x} - \underline{x}$
  • 中点（midpoint）：$\text{mid}(\mathbf{x}) = \dfrac{\underline{x} + \overline{x}}{2}$
  • 半径（radius）：$\text{rad}(\mathbf{x}) = \dfrac{w(\mathbf{x})}{2} = \dfrac{\overline{x} - \underline{x}}{2}$
  • 絶対値（absolute value）：$|\mathbf{x}| = \max(|\underline{x}|, |\overline{x}|)$
  • 最小絶対値（mignitude）：$\langle\mathbf{x}\rangle = \begin{cases} \min(|\underline{x}|, |\overline{x}|) & \text{if } 0 \notin \mathbf{x} \\ 0 & \text{otherwise} \end{cases}$

特に $w(\mathbf{x}) = 0$ のとき、$\mathbf{x}$ は点区間（point interval、thin interval）であり、 通常の実数と同一視できる。

四則演算

区間の四則演算は、演算結果として取りうる値の集合全体を包含する最小の区間として定義される。 形式的には、二項演算 $\circ \in \{+, -, \times, /\}$ に対し：

$$\mathbf{x} \circ \mathbf{y} = \{ x \circ y \mid x \in \mathbf{x}, \, y \in \mathbf{y} \}$$

の最小包含区間（interval hull）として定義される。以下、各演算の具体的な公式を与える。

定義（区間の加法）

$$[\underline{x}, \overline{x}] + [\underline{y}, \overline{y}] = [\underline{x} + \underline{y},\; \overline{x} + \overline{y}]$$

証明

任意の $a \in [\underline{x}, \overline{x}]$ と $b \in [\underline{y}, \overline{y}]$ をとる。 定義より $\underline{x} \leq a \leq \overline{x}$ かつ $\underline{y} \leq b \leq \overline{y}$ であるから、 辺々加えて：

$$\underline{x} + \underline{y} \leq a + b \leq \overline{x} + \overline{y}$$

したがって $a + b \in [\underline{x} + \underline{y}, \overline{x} + \overline{y}]$ である。 逆に、$a = \underline{x}, b = \underline{y}$ とすれば下端が達成され、 $a = \overline{x}, b = \overline{y}$ とすれば上端が達成されるので、 これが最小包含区間である。$\square$

定義（区間の減法）

$$[\underline{x}, \overline{x}] - [\underline{y}, \overline{y}] = [\underline{x} - \overline{y},\; \overline{x} - \underline{y}]$$

証明

$a \in [\underline{x}, \overline{x}]$, $b \in [\underline{y}, \overline{y}]$ のとき、 $-\overline{y} \leq -b \leq -\underline{y}$ であるから：

$$\underline{x} - \overline{y} \leq a - b \leq \overline{x} - \underline{y}$$

同様に端点で等号が成立し、最小包含区間である。$\square$

定義（区間の乗法）

$S = \{\underline{x}\,\underline{y},\; \underline{x}\,\overline{y},\; \overline{x}\,\underline{y},\; \overline{x}\,\overline{y}\}$ として：

$$[\underline{x}, \overline{x}] \times [\underline{y}, \overline{y}] = [\min S,\; \max S]$$

証明

$f(a, b) = ab$ は各変数について線形であるから、 矩形領域 $[\underline{x}, \overline{x}] \times [\underline{y}, \overline{y}]$ 上での最大・最小は 必ず頂点（端点の組み合わせ）で達成される。 頂点は $(\underline{x}, \underline{y})$, $(\underline{x}, \overline{y})$, $(\overline{x}, \underline{y})$, $(\overline{x}, \overline{y})$ の4つであるから、$S$ の最小値と最大値がそれぞれ下端と上端を与える。$\square$

乗法の特殊ケース

両区間の符号が確定している場合、4つの積をすべて計算する必要はない：

  • $\underline{x} \geq 0$ かつ $\underline{y} \geq 0$（両方非負）のとき： $\mathbf{x} \times \mathbf{y} = [\underline{x}\,\underline{y},\; \overline{x}\,\overline{y}]$
  • $\overline{x} \leq 0$ かつ $\overline{y} \leq 0$（両方非正）のとき： $\mathbf{x} \times \mathbf{y} = [\overline{x}\,\overline{y},\; \underline{x}\,\underline{y}]$
  • $\underline{x} \geq 0$ かつ $\overline{y} \leq 0$（$\mathbf{x}$ 非負、$\mathbf{y}$ 非正）のとき： $\mathbf{x} \times \mathbf{y} = [\overline{x}\,\underline{y},\; \underline{x}\,\overline{y}]$

実装では符号による場合分けを行い、不要な乗算を省略することで高速化を図る。

定義（区間の除法）

$0 \notin [\underline{y}, \overline{y}]$ のとき：

$$[\underline{x}, \overline{x}] \;/\; [\underline{y}, \overline{y}] = [\underline{x}, \overline{x}] \times \left[\dfrac{1}{\overline{y}},\; \dfrac{1}{\underline{y}}\right]$$

証明

$0 \notin [\underline{y}, \overline{y}]$ のとき、$b \in [\underline{y}, \overline{y}]$ なら $b \neq 0$ である。 $g(b) = 1/b$ は $0$ を含まない区間上で単調であるから：

  • $\underline{y} > 0$ のとき：$g$ は単調減少、$1/b \in [1/\overline{y}, 1/\underline{y}]$
  • $\overline{y} < 0$ のとき：$g$ は単調減少、$1/b \in [1/\overline{y}, 1/\underline{y}]$

いずれの場合も $\{1/b \mid b \in [\underline{y}, \overline{y}]\} = [1/\overline{y}, 1/\underline{y}]$ となる。 よって除法は逆数の区間との乗法に帰着される。$\square$

ゼロ除算について

$0 \in [\underline{y}, \overline{y}]$ の場合、逆数 $1/b$ は $b \to 0$ で発散するため、 通常の区間演算では除法は未定義となる。 拡張区間演算（extended interval arithmetic）では、 $[-\infty, +\infty]$ を含む拡張区間を導入してこのケースを扱う。 例えば $[1,2] / [-1, 1] = [-\infty, +\infty]$ のように、 結果として全実数を返す。より精密には、 $\underline{y} = 0$ かつ $\overline{y} > 0$ のとき $[1,2] / [0, 1] = [1, +\infty]$ のように片側無限区間を返すこともある。

具体的な数値例

例：加法

$$[1.5, 2.5] + [3.0, 4.0] = [1.5 + 3.0,\; 2.5 + 4.0] = [4.5, 6.5]$$

例：乗法

$\mathbf{x} = [1.5, 2.5]$, $\mathbf{y} = [3.0, 4.0]$ とする。 両区間とも正であるから：

$$[1.5, 2.5] \times [3.0, 4.0] = [1.5 \times 3.0,\; 2.5 \times 4.0] = [4.5, 10.0]$$

例：符号が混在する乗法

$\mathbf{x} = [-1, 2]$, $\mathbf{y} = [-3, 4]$ とする。端点の4つの積は：

$$S = \{(-1)(-3),\; (-1)(4),\; (2)(-3),\; (2)(4)\} = \{3, -4, -6, 8\}$$ $$[-1, 2] \times [-3, 4] = [\min S, \max S] = [-6, 8]$$

区間演算の包含性定理

定理（区間演算の基本定理：包含性）

$f(x_1, \ldots, x_n)$ を変数 $x_1, \ldots, x_n$ の有理式（四則演算のみで構成される式）とする。 $\mathbf{f}(\mathbf{x}_1, \ldots, \mathbf{x}_n)$ を、$f$ の各演算を対応する区間演算に置き換えた 区間拡張（interval extension）とする。

各 $x_i \in \mathbf{x}_i$ ($i = 1, \ldots, n$) ならば（除法で $0$ 除算が生じない限り）：

$$f(x_1, \ldots, x_n) \in \mathbf{f}(\mathbf{x}_1, \ldots, \mathbf{x}_n)$$

証明

$f$ の構造に関する帰納法で示す。

基底：$f(x_i) = x_i$ のとき、$x_i \in \mathbf{x}_i$ は自明。 $f = c$（定数）のとき、$c \in [c, c]$ も自明。

帰納ステップ：$f = g \circ h$（$\circ \in \{+,-,\times,/\}$）のとき、 帰納法の仮定より $g(x_1, \ldots, x_n) \in \mathbf{g}(\mathbf{x}_1, \ldots, \mathbf{x}_n)$ かつ $h(x_1, \ldots, x_n) \in \mathbf{h}(\mathbf{x}_1, \ldots, \mathbf{x}_n)$。 区間演算の定義より、$a \in \mathbf{a}$, $b \in \mathbf{b}$ ならば $a \circ b \in \mathbf{a} \circ \mathbf{b}$ であるから：

$$f(x_1, \ldots, x_n) = g(\ldots) \circ h(\ldots) \in \mathbf{g}(\ldots) \circ \mathbf{h}(\ldots) = \mathbf{f}(\mathbf{x}_1, \ldots, \mathbf{x}_n)$$

$\square$

注意：包含性の意味

包含性は「真の値が結果区間に必ず含まれる」ことを保証するが、 結果区間が最小（最も狭い）であることは一般には保証しない。 同じ変数が複数回出現する場合、各出現が独立に扱われるため、 結果区間が真の値域よりも大きくなることがある。 これが次節で扱う「依存性問題」である。

依存性問題

区間演算は包含性を保証する一方で、過大評価（overestimation）という本質的な問題を抱えている。 これは依存性問題（dependency problem）あるいはwrapping effectと呼ばれ、 同一変数の複数回の出現を区間演算が独立に扱ってしまうことに起因する。

具体例

例1：$x - x$ の過大評価

$x \in [0, 1]$ のとき、実際には $x - x = 0$ が常に成り立つ。しかし区間演算では：

$$[0, 1] - [0, 1] = [0 - 1,\; 1 - 0] = [-1, 1]$$

幅2の区間が得られ、真の値域 $\{0\}$ を大きく過大評価している。 これは、減法 $\mathbf{x} - \mathbf{y}$ において $\mathbf{x}$ と $\mathbf{y}$ が 「同じ変数である」という情報が失われ、 $\mathbf{x}$ の最大値から $\mathbf{y}$ の最小値を引く組み合わせ（$1 - 0 = 1$）が許容されてしまうためである。

例2：式の形による結果の違い

$x \in [2, 3]$ のとき、$f(x) = x^2 - 2x + 1 = (x-1)^2$ を考える。真の値域は $[1, 4]$ である。

直接計算（$x$ が3回出現）：

$$\begin{align} [2,3]^2 - 2 \cdot [2,3] + 1 &= [4, 9] - [4, 6] + [1, 1] \\ &= [4-6,\; 9-4] + [1, 1] \\ &= [-2, 5] + [1, 1] \\ &= [-1, 6] \end{align}$$

真の値域 $[1, 4]$ に対し $[-1, 6]$ は幅7の過大評価となっている。

変形後の計算（$x$ が1回のみ出現）：

$$([2,3] - [1,1])^2 = [1, 2]^2 = [1, 4]$$

真の値域 $[1, 4]$ と一致する。$x$ の出現回数が1回に減ったため、依存性問題が完全に解消された。

根本原因

依存性問題の根本原因は、区間演算が各変数の各出現を独立な変数として扱うことにある。 数学的に言えば、区間演算は集合演算

$$\mathbf{x} \circ \mathbf{y} = \{ x \circ y \mid x \in \mathbf{x}, \, y \in \mathbf{y} \}$$

を忠実に計算するが、「$x$ と $y$ は実は同じ変数の同じ値である」という制約を取り込む仕組みがない。 したがって、式中に同じ変数が多く出現するほど過大評価が大きくなる傾向がある。

対策

1. 式の書き換え

上の例2で示したように、数学的に等価な式に変形して変数の出現回数を減らすことで過大評価を軽減できる。 特に、$f(x) = x^2 - 2x + 1$ を $f(x) = (x-1)^2$ と変形すれば $x$ は1回しか出現しないため、 依存性問題は生じない。 ただし、すべての式に対してこのような変形が可能とは限らない。

2. アフィン演算（Affine Arithmetic）

定義（アフィン形式）

アフィン形式（affine form）とは、以下の形の式である：

$$\hat{x} = x_0 + x_1 \varepsilon_1 + x_2 \varepsilon_2 + \cdots + x_n \varepsilon_n$$

ここで $x_0 \in \mathbb{R}$ は中心値、$x_i \in \mathbb{R}$ ($i \geq 1$) は偏差係数、 $\varepsilon_i \in [-1, 1]$ はノイズシンボル（noise symbol）と呼ばれる独立な変数である。

アフィン演算の核心は、異なる量に含まれる同じノイズシンボルが依存関係を表現する点にある。 例えば、$\hat{x} = 5 + 2\varepsilon_1$ かつ $\hat{y} = \hat{x}$ とすれば、 $\hat{x} - \hat{y} = (5 + 2\varepsilon_1) - (5 + 2\varepsilon_1) = 0$ となり、 通常の区間演算では過大評価される $x - x$ が正しく $0$ と計算される。

アフィン演算の四則演算

$\hat{x} = x_0 + \displaystyle\sum_i x_i \varepsilon_i$, $\hat{y} = y_0 + \displaystyle\sum_i y_i \varepsilon_i$ に対する演算規則は以下の通りである。

演算 結果のアフィン形式 精度
$\hat{x} + \hat{y}$ $(x_0 + y_0) + \displaystyle\sum_i (x_i + y_i)\varepsilon_i$ 厳密（新シンボル不要）
$\alpha \hat{x} + \beta$ $(\alpha x_0 + \beta) + \displaystyle\sum_i \alpha x_i \varepsilon_i$ 厳密
$\hat{x} \cdot \hat{y}$ $x_0 y_0 + \displaystyle\sum_i (x_0 y_i + y_0 x_i)\varepsilon_i + z_{n+1}\varepsilon_{n+1}$ 近似（2次項を新シンボルに集約）

乗算では2次の項 $\displaystyle\sum_{i,j} x_i y_j \varepsilon_i \varepsilon_j$ がアフィン形式で表現できないため、 その絶対値の上界 $z_{n+1} = \displaystyle\sum_i |x_i| \cdot \displaystyle\sum_j |y_j|$ として 新しいノイズシンボル $\varepsilon_{n+1}$ に吸収する。

アフィン演算 vs 区間演算の比較

$x \in [2, 4]$ に対し $f(x) = x^2 - x$ を計算する。

  • 区間演算：$[2,4]^2 - [2,4] = [4,16] - [2,4] = [0, 14]$（真の値域 $[2, 12]$ に対し過大評価）
  • アフィン演算：$\hat{x} = 3 + \varepsilon_1$ とおくと、 $\hat{x}^2 \approx 9 + 6\varepsilon_1 + \varepsilon_2$（$|\varepsilon_2| \leq 1$）、 $\hat{x}^2 - \hat{x} = 6 + 5\varepsilon_1 + \varepsilon_2 \in [0, 12]$ で、より正確

アフィン演算の限界と改良

アフィン演算は線形演算に対しては依存性を正確に追跡できるが、 非線形演算では新しいノイズシンボルが導入されるため、長い計算ではシンボル数が増大する。 改良版として修正アフィン演算（Modified Affine Arithmetic, MAA）や 一般化アフィン演算（GAA）が提案されており、 非線形演算でのシンボル管理を改善する。

3. Taylor模型（Taylor Models）

Taylor模型は、関数の値をTaylor多項式と区間剰余項の組で表現する手法である：

$$f(x) \in T_k(x) + \mathbf{r}$$

ここで $T_k(x)$ は $k$ 次Taylor多項式、$\mathbf{r}$ は剰余項を包含する区間である。 多項式部分で関数の振る舞いを精密に捉え、区間部分は剰余のみを包含すればよいため、 通常の区間演算と比べて過大評価が大幅に抑制される。 Makino & Berz (1996) により体系化され、ビーム物理学や天体力学で応用されている。

Taylor模型の演算

Taylor模型 $\mathcal{T}_f = (T_f, \mathbf{r}_f)$ と $\mathcal{T}_g = (T_g, \mathbf{r}_g)$ に対する演算は以下の通りである。

  • 加算：$\mathcal{T}_f + \mathcal{T}_g = (T_f + T_g,\; \mathbf{r}_f + \mathbf{r}_g)$
  • 乗算：$T_f \cdot T_g$ の $k$ 次までの項を多項式部分に、$k+1$ 次以上の項の範囲を区間剰余に加算
  • 合成：$f(g(x))$ の Taylor 展開を $k$ 次まで自動微分で計算し、剰余を区間で包含
手法 過大評価の次数 計算コスト 適用場面
区間演算 $O(w(\mathbf{x}))$ 最小 狭い区間での基本演算
アフィン演算 $O(w(\mathbf{x})^2)$ 中程度 線形性の高い計算
$k$ 次 Taylor 模型 $O(w(\mathbf{x})^{k+1})$ $O(n^k)$ ODE の長時間積分、大域最適化

Taylor模型の応用例

  • 常微分方程式の精度保証：Taylor 模型をステップごとに更新し、解の存在と一意性を保証しつつ長時間積分する（COSY INFINITY, CAPD）
  • 大域最適化：Taylor 模型の多項式部分から関数の下界を効率的に計算し、分枝限定法を加速する
  • 宇宙力学：小惑星の軌道不確定性を Taylor 模型で伝播させ、衝突確率を厳密に評価する

4. 平均値形式（Mean Value Form）

平均値形式

$f$ が微分可能な関数、$m = \text{mid}(\mathbf{x})$ のとき：

$$f(\mathbf{x}) \subseteq f(m) + f'(\mathbf{x}) \cdot (\mathbf{x} - m)$$

右辺は通常の区間演算による包含で計算される。 $\mathbf{x}$ の幅が小さいとき、$(\mathbf{x} - m)$ の幅は $w(\mathbf{x})/2$ 程度であるため、 直接の区間評価よりも狭い結果が得られることが多い。

根拠

平均値の定理より、任意の $x \in \mathbf{x}$ に対し、ある $\xi \in \mathbf{x}$ が存在して $f(x) = f(m) + f'(\xi)(x - m)$ が成り立つ。 $\xi \in \mathbf{x}$ であるから $f'(\xi) \in f'(\mathbf{x})$ であり、 $x - m \in \mathbf{x} - m$ であるから：

$$f(x) = f(m) + f'(\xi)(x - m) \in f(m) + f'(\mathbf{x}) \cdot (\mathbf{x} - m)$$

$\square$

区間ニュートン法

通常のニュートン法の復習

$f(x^*) = 0$ の根 $x^*$ を求めるためのニュートン法（Newton's method）は、 初期値 $x_0$ から始めて以下の反復を行う：

$$x_{n+1} = x_n - \dfrac{f(x_n)}{f'(x_n)}$$

通常のニュートン法は収束が速い（二次収束）が、以下の限界がある：

  • 初期値が根の十分近くにないと収束しない可能性がある
  • 収束したとしても、他の根を見逃す可能性がある
  • 根の存在を数学的に証明する機能がない

区間ニュートン作用素

定義（区間ニュートン作用素）

$f$ が連続微分可能、$\mathbf{x} = [\underline{x}, \overline{x}]$、 $m = \text{mid}(\mathbf{x})$ のとき、区間ニュートン作用素を以下で定義する：

$$N(\mathbf{x}) = m - \dfrac{f(m)}{f'(\mathbf{x})}$$

ここで $f'(\mathbf{x})$ は $\mathbf{x}$ 上の $f'$ の区間評価であり、除法は区間除法である。

反復は、この作用素と元の区間の交差をとることで行われる：

$$\mathbf{x}_{\text{new}} = \mathbf{x} \cap N(\mathbf{x})$$

三つの結果

区間ニュートン作用素の結果には、以下の3つの場合がある。 これらが区間ニュートン法の真の威力である。

定理（区間ニュートン法の判定）

$f$ は $\mathbf{x}$ 上で連続微分可能、$0 \notin f'(\mathbf{x})$ とする。このとき：

  1. $N(\mathbf{x}) \cap \mathbf{x} = \emptyset$ ならば：$\mathbf{x}$ 内に $f$ の根は存在しない（不在証明）
  2. $N(\mathbf{x}) \subseteq \text{int}(\mathbf{x})$（真の包含）ならば：$\mathbf{x}$ 内に $f$ の根がちょうど一つ存在する（存在・一意性の証明）
  3. それ以外：根が存在するかどうかは不明、区間が絞り込まれる

証明のスケッチ

ケース1（不在証明）：

$f(x^*) = 0$ かつ $x^* \in \mathbf{x}$ と仮定する。 平均値の定理より、ある $\xi \in \mathbf{x}$ が存在して：

$$f(m) = f(m) - f(x^*) = f'(\xi)(m - x^*)$$ $$x^* = m - \dfrac{f(m)}{f'(\xi)}$$

$\xi \in \mathbf{x}$ であるから $f'(\xi) \in f'(\mathbf{x})$、 したがって $x^* \in m - f(m)/f'(\mathbf{x}) = N(\mathbf{x})$ である。 同時に $x^* \in \mathbf{x}$ であるから $x^* \in N(\mathbf{x}) \cap \mathbf{x}$。 対偶をとれば、$N(\mathbf{x}) \cap \mathbf{x} = \emptyset$ ならば $\mathbf{x}$ に根は存在しない。

ケース2（存在・一意性の証明）：

$N(\mathbf{x}) \subseteq \text{int}(\mathbf{x})$ のとき、 写像 $g(x) = x - f(x)/f'(c)$（$c$ は $f'(\mathbf{x})$ の適切な代表元）が $\mathbf{x}$ から $\mathbf{x}$ への縮小写像となることを示せる。 Brouwer の不動点定理（あるいは縮小写像の原理）より、 $\mathbf{x}$ 内にちょうど一つの不動点、すなわち $f$ の根が存在する。$\square$

区間ニュートン法の威力

通常のニュートン法は根に収束する「数値的な手法」に過ぎないが、 区間ニュートン法は根の存在を証明し、不在を証明できる 数学的証明の道具である。 これは数値解析において極めて稀な性質であり、 計算機援用証明の中核技術となっている。

具体例

例：$f(x) = x^2 - 2$ の根を $[1, 2]$ で求める

$f(x) = x^2 - 2$、$f'(x) = 2x$ である。初期区間 $\mathbf{x} = [1, 2]$ とする。

第1反復：

  • $m = \text{mid}([1, 2]) = 1.5$
  • $f(m) = f(1.5) = 2.25 - 2 = 0.25$
  • $f'(\mathbf{x}) = 2 \cdot [1, 2] = [2, 4]$
  • $N([1, 2]) = 1.5 - \dfrac{0.25}{[2, 4]} = 1.5 - [0.0625, 0.125] = [1.375, 1.4375]$

$[1.375, 1.4375] \subset [1, 2]$（真の包含）であるから、定理のケース2より、$[1, 2]$ 内に $f$ の根がちょうど一つ存在することが証明された。

実際、$\sqrt{2} \approx 1.41421356\ldots$ はこの区間 $[1.375, 1.4375]$ に含まれる。

第2反復：$\mathbf{x} = [1.375, 1.4375]$ として反復を続ける。

  • $m = 1.40625$
  • $f(m) = 1.40625^2 - 2 = 1.97754 - 2 = -0.02246$
  • $f'(\mathbf{x}) = 2 \cdot [1.375, 1.4375] = [2.75, 2.875]$
  • $N(\mathbf{x}) = 1.40625 - \dfrac{-0.02246}{[2.75, 2.875]} = 1.40625 + [0.00781, 0.00817] = [1.41406, 1.41442]$

幅が $0.0625$ から $0.00036$ へと大幅に縮小し、$\sqrt{2}$ に急速に収束している。

二分法との組み合わせ

実用的には、区間ニュートン法は二分法（bisection）と組み合わせて使われる。 探索区間を二分割し、各部分区間に対して区間ニュートン法を適用する：

  1. 初期区間 $\mathbf{x}$ を二分して $\mathbf{x}_L$ と $\mathbf{x}_R$ を得る
  2. 各部分区間に対し区間ニュートン作用素を計算する
  3. $N(\mathbf{x}_i) \cap \mathbf{x}_i = \emptyset$ ならその部分区間を棄却する（根なしが証明された）
  4. $N(\mathbf{x}_i) \subseteq \text{int}(\mathbf{x}_i)$ ならその部分区間に根が1つ存在することが確定
  5. それ以外の場合、さらに二分して再帰的に処理する

この方法により、ある区間内のすべての根を漏れなく求めることができ、 しかも各根の存在が数学的に証明される。 通常のニュートン法では得られない、大域的かつ厳密な根の探索が可能となる。

初等関数の区間拡張

四則演算に加え、初等関数（冪乗、指数・対数、三角関数など）を区間に対して適用する方法が必要である。 基本原理は同じで、関数の値域の最小包含区間を求める。 単調関数に対しては端点評価のみで済むが、 非単調関数（三角関数など）では区間内での極値を考慮する必要がある。

単調関数の場合

定理（単調関数の区間拡張）

$f$ が区間 $[\underline{x}, \overline{x}]$ 上で単調増加であるとき：

$$f([\underline{x}, \overline{x}]) = [f(\underline{x}), f(\overline{x})]$$

$f$ が単調減少であるとき：

$$f([\underline{x}, \overline{x}]) = [f(\overline{x}), f(\underline{x})]$$

この定理は中間値の定理と単調性から直ちに従う。 指数関数 $\exp$、対数関数 $\log$、双曲線関数 $\sinh$ などは全域で単調であるため、端点評価のみで最小包含区間が得られる。

例：指数関数の区間拡張

$\exp$ は単調増加であるから：

$$\exp([1, 2]) = [\exp(1), \exp(2)] = [2.71828\ldots, 7.38905\ldots]$$

浮動小数点で実装する場合は、下端を $-\infty$ 方向に丸め、上端を $+\infty$ 方向に丸める。

冪乗 $x^n$ の場合

区間の冪乗

整数冪 $[\underline{x}, \overline{x}]^n$ の計算は $n$ の偶奇と区間の符号に依存する：

  • $n$ が奇数のとき：$x^n$ は単調増加であるから $[\underline{x}^n, \overline{x}^n]$
  • $n$ が偶数のとき：$x^n$ は $x = 0$ で最小値をとり、以下の場合分けが必要：
    • $\underline{x} \geq 0$（全体が非負）：$[\underline{x}^n, \overline{x}^n]$
    • $\overline{x} \leq 0$（全体が非正）：$[\overline{x}^n, \underline{x}^n]$
    • $\underline{x} < 0 < \overline{x}$（ゼロを跨ぐ）：$[0, \max(\underline{x}^n, \overline{x}^n)]$

例：$[-1, 2]^2$ と $[-1, 2] \times [-1, 2]$ の違い

冪乗としての $x^2$（$x$ が同じ値）：

$$[-1, 2]^2 = [0, \max((-1)^2, 2^2)] = [0, 4]$$

乗法としての $\mathbf{x} \times \mathbf{x}$（依存性を無視）：

$$[-1, 2] \times [-1, 2] = [\min\{1, -2, -2, 4\}, \max\{1, -2, -2, 4\}] = [-2, 4]$$

冪乗の方が真の値域 $[0, 4]$ と一致する。これは $x \times x$ の計算で $x$ の2回の出現を独立に扱ってしまう依存性問題の典型例である。 区間演算ライブラリでは $x^2$ を乗法ではなく専用の冪乗関数として実装することで、この過大評価を回避する。

三角関数の場合

三角関数は非単調であるため、区間拡張にはより複雑な処理が必要である。 $\sin$ を例に説明する。

$\sin$ の区間拡張

$\sin([\underline{x}, \overline{x}])$ を計算するには、区間 $[\underline{x}, \overline{x}]$ 内に $\sin$ の極値点（$x = \pi/2 + k\pi$, $k \in \mathbb{Z}$）が含まれるかを判定する必要がある：

  1. $[\underline{x}, \overline{x}]$ が $\sin$ の極大点（$\sin = 1$）を含むかチェック
  2. $[\underline{x}, \overline{x}]$ が $\sin$ の極小点（$\sin = -1$）を含むかチェック
  3. 極値点を含まない場合は端点値のみで判定可能（$\sin$ は部分区間上で単調）
  4. 上端は $\max(\sin\underline{x}, \sin\overline{x})$（極大点を含む場合は $1$）
  5. 下端は $\min(\sin\underline{x}, \sin\overline{x})$（極小点を含む場合は $-1$）

例：$\sin([0, \pi])$

区間 $[0, \pi]$ 内に極大点 $x = \pi/2$ が含まれる。端点での値は $\sin(0) = 0$, $\sin(\pi) = 0$。 極大点での値は $\sin(\pi/2) = 1$。よって：

$$\sin([0, \pi]) = [\min(0, 0), 1] = [0, 1]$$

もしこれを端点のみで計算すると $[\min(0,0), \max(0,0)] = [0, 0]$ と誤った結果になる。 区間内の極値点の検出が不可欠であることが分かる。

実装上の注意

三角関数の区間拡張では、極値点の位置判定に $\pi$ の値が必要であるが、 $\pi$ は無理数であるため厳密には表現できない。 実装では $\pi$ の区間包含 $[\underline{\pi}, \overline{\pi}]$ を用いて安全に判定を行う必要がある。 この「$\pi$ の問題」は区間演算ライブラリの実装上の重要な技術的課題の一つである。

区間行列演算

区間演算のスカラーからベクトル・行列への拡張は、連立方程式の精度保証（第11章）をはじめとする多くの応用の基盤となる。 区間行列とは、各成分が区間である行列のことである。

定義（区間ベクトルと区間行列）

区間ベクトル $[\mathbf{x}] = ([\underline{x}_i, \overline{x}_i])_{i=1}^n \in \mathbb{IR}^n$ は、 各成分が区間であるベクトルである。幾何学的には $\mathbb{R}^n$ 内の軸平行な直方体（box）を表す。

区間行列 $[\mathbf{A}] = ([\underline{a}_{ij}, \overline{a}_{ij}])_{i,j} \in \mathbb{IR}^{m \times n}$ は、各成分が区間である行列である。

区間行列の演算

区間行列の加法・乗法は、通常の行列演算をスカラー演算の組み合わせに分解し、各スカラー演算を区間演算に置き換えたものである。

区間行列の乗法

$[\mathbf{A}] \in \mathbb{IR}^{m \times p}$, $[\mathbf{B}] \in \mathbb{IR}^{p \times n}$ に対して：

$$([\mathbf{A}][\mathbf{B}])_{ij} = \displaystyle\sum_{k=1}^{p} [\mathbf{A}]_{ik} \cdot [\mathbf{B}]_{kj}$$

各成分は $p$ 個の区間積の区間和として計算される。包含性はスカラーの場合から自動的に従う。

過大評価の問題

区間行列の乗法では、スカラーの場合と同様に依存性問題による過大評価が生じる。 特に内積 $\displaystyle\sum_{k} a_{ik} b_{kj}$ の計算では、各項の和を逐次的に取るため過大評価が蓄積する。

$n \times n$ の区間行列の乗法のコストは $O(n^3)$ 回の区間乗算と区間加算であり、 通常の行列乗算の約4倍のコスト（各区間演算が2つの浮動小数点演算を要するため）である。

区間連立一次方程式

区間行列 $[\mathbf{A}] \in \mathbb{IR}^{n \times n}$ と区間ベクトル $[\mathbf{b}] \in \mathbb{IR}^n$ に対して、 解集合を考える：

$$\Sigma([\mathbf{A}], [\mathbf{b}]) = \{ \mathbf{x} \in \mathbb{R}^n \mid A\mathbf{x} = \mathbf{b} \text{ for some } A \in [\mathbf{A}], \, \mathbf{b} \in [\mathbf{b}] \}$$

解集合の構造

一般に解集合 $\Sigma$ は凸集合ではないため、区間ベクトル（直方体）で包含すると過大評価が生じる。 Oettli-Prager (1964) の定理は、解集合の包含に関する以下の特徴づけを与える：

定理（Oettli-Prager, 1964）

$A_c = \text{mid}([\mathbf{A}])$, $\Delta = \text{rad}([\mathbf{A}])$, $b_c = \text{mid}([\mathbf{b}])$, $\delta = \text{rad}([\mathbf{b}])$ とする。 このとき $\mathbf{x} \in \Sigma([\mathbf{A}], [\mathbf{b}])$ であるための必要十分条件は：

$$|A_c \mathbf{x} - b_c| \leq \Delta |\mathbf{x}| + \delta$$

（不等式は成分ごとに成立）。ここで $|A_c \mathbf{x} - b_c|$ は残差の絶対値ベクトルである。

この定理は、解集合が行列とベクトルの不確かさの範囲内での残差条件によって特徴づけられることを示している。 区間連立一次方程式の直接的な求解（解集合の最小包含区間の計算）はNP困難であることが知られており、 実用的には第11章で述べるKrawczyk法やRump法が用いられる。

制約伝播

制約伝播（constraint propagation）は、区間演算と組み合わせて用いられる手法であり、 等式・不等式制約から区間の縮小（絞り込み）を行う。 大域的最適化、非線形方程式の求解、ロバスト制御設計などで重要な役割を果たす。

基本的なアイデア

変数 $x, y, z$ に対する制約 $z = x + y$ が与えられ、 各変数の区間が $x \in [1, 5]$, $y \in [2, 6]$, $z \in [4, 8]$ であるとする。 区間演算による前進的な評価では $z \in [1+2, 5+6] = [3, 11]$ が得られるが、 $z$ の既知の区間 $[4, 8]$ と交差をとれば $z \in [4, 8]$ に絞られる。

さらに、$z \in [4, 8]$ と $y \in [2, 6]$ から逆方向に $x$ の区間を絞ることもできる。 $x = z - y$ より $x \in [4-6, 8-2] = [-2, 6]$。既知の $x \in [1, 5]$ との交差で $x \in [1, 5]$（この場合は縮小なし）。 同様に $y = z - x$ より $y \in [4-5, 8-1] = [-1, 7] \cap [2, 6] = [2, 6]$。

制約伝播のアルゴリズム

制約の集合 $\{c_1, c_2, \ldots, c_m\}$ と変数の初期区間 $[\mathbf{x}_1], \ldots, [\mathbf{x}_n]$ が与えられたとき：

  1. キューにすべての制約を入れる
  2. キューから制約 $c_k$ を取り出し、関連する変数の区間を前進・後退の両方向で縮小する
  3. いずれかの変数の区間が縮小されたら、その変数に関連する他の制約をキューに追加する
  4. キューが空になるか、縮小量が閾値以下になれば終了する
  5. いずれかの変数の区間が空になれば、制約系は実行不可能（解なし）と判定できる

基本演算に対する縮小規則

四則演算の逆方向伝播

制約 $z = x \circ y$ に対する各変数の縮小規則：

制約 $z$ の縮小 $x$ の縮小 $y$ の縮小
$z = x + y$ $[z] \cap ([x] + [y])$ $[x] \cap ([z] - [y])$ $[y] \cap ([z] - [x])$
$z = x \times y$ $[z] \cap ([x] \times [y])$ $[x] \cap ([z] / [y])$ $[y] \cap ([z] / [x])$
$z = x^2$ $[z] \cap [x]^2$ $[x] \cap ([\pm\sqrt{\underline{z}}, \pm\sqrt{\overline{z}}])$（符号考慮）

乗法の逆方向伝播で除法を用いる場合、$0 \in [y]$ のときは拡張区間除法が必要となることに注意する。

応用例：大域的最適化の加速

例：制約伝播による探索空間の縮小

問題：$f(x, y) = x^2 + y^2$ を $x + y \geq 3$, $x \in [0, 5]$, $y \in [0, 5]$ のもとで最小化する。

制約伝播：制約 $x + y \geq 3$ から $x \geq 3 - y$ より $x \geq 3 - 5 = -2$（既知の $x \geq 0$ と合わせて変化なし）。 しかし $y \geq 3 - x$ より $y \geq 3 - 5 = -2$（同様に変化なし）。

関数値の上界 $f^*_{\text{UB}}$ が例えば $10$ と判明しているとする。 すると $x^2 \leq 10$ より $x \in [0, \sqrt{10}] \cap [0, 5] = [0, 3.162]$ に縮小される。 同様に $y \in [0, 3.162]$。

さらに制約 $x + y \geq 3$ と $y \leq 3.162$ から $x \geq 3 - 3.162 = -0.162$（変化なし）。 $x \leq 3.162$ から $y \geq 3 - 3.162 = -0.162$（変化なし）。

探索空間が $[0,5]^2$（面積25）から $[0, 3.162]^2$（面積10）に縮小された。区間分枝限定法と組み合わせることで、探索の効率が大幅に向上する。

HC4 アルゴリズム

HC4（Benhamou et al., 1999）は、制約伝播の実用的なアルゴリズムである。 制約式の構文木（DAG）に沿って前進評価と後退縮小を交互に行う。 各制約について $O(k)$（$k$ は制約式のサイズ）のコストで1回の伝播が可能であり、 不動点に達するまで反復する。制約充足問題（CSP）のソルバーの中核として広く実装されている。

向き付き丸め

実装上の課題

区間演算の包含性を計算機上で厳密に保証するためには、 区間の下端と上端を浮動小数点数で表現する際の丸めの方向を制御する必要がある。 例えば、$[1, 1] / [3, 3]$ を計算すると真の結果は $[1/3, 1/3]$ であるが、 $1/3$ は有限桁の浮動小数点数では正確に表現できない。 通常の「最近接丸め」を使うと、下端と上端が同じ浮動小数点数に丸められ、 真の値 $1/3$ がその区間に含まれない恐れがある。

定義（向き付き丸め）

区間演算の実装では、以下の向き付き丸め（directed rounding）を用いる：

  • 下端の計算：$-\infty$ 方向に丸める（round toward $-\infty$, floor） $$\underline{z} = \text{RD}(f(\underline{x}, \overline{x}, \underline{y}, \overline{y}))$$
  • 上端の計算：$+\infty$ 方向に丸める（round toward $+\infty$, ceiling） $$\overline{z} = \text{RU}(g(\underline{x}, \overline{x}, \underline{y}, \overline{y}))$$

ここで $\text{RD}$ は下方丸め、$\text{RU}$ は上方丸めを表す。 これにより、浮動小数点で表現された区間 $[\underline{z}, \overline{z}]$ が 真の結果区間を必ず包含することが保証される。

IEEE 754 と丸めモード

IEEE 754 浮動小数点規格は、以下の4つの丸めモードを定義している：

  1. 最近接偶数丸め（round to nearest, ties to even）：既定の丸めモード
  2. $+\infty$ 方向丸め（round toward $+\infty$）：区間上端に使用
  3. $-\infty$ 方向丸め（round toward $-\infty$）：区間下端に使用
  4. ゼロ方向丸め（round toward zero）：切り捨て

区間演算で必要なのは主に2番と3番である。

ハードウェアサポート

ほとんどの現代的な CPU は、FPU（浮動小数点演算ユニット）の制御ワードを変更することで 丸めモードを切り替える機能を備えている。 例えば x86 アーキテクチャでは、SSE の MXCSR レジスタにより丸めモードを制御できる。

ただし、丸めモードの切り替えにはオーバーヘッドがあり、 頻繁な切り替えはパイプラインのストールを引き起こす可能性がある。 このため、実用的な区間演算ライブラリでは以下のような最適化が行われる：

  • 下端の計算をまとめて行い（$-\infty$ 丸めモード）、次に上端をまとめて行う（$+\infty$ 丸めモード）
  • あるいは、常に上方丸めモードで計算し、下端は符号を反転して計算する （$\text{RD}(a + b) = -\text{RU}(-a - b)$ を利用）

ソフトウェア的代替手法

ハードウェアの丸めモード変更が利用できない場合や、 より移植性の高い実装が必要な場合には、ソフトウェア的な手法が用いられる：

  • 1 ULP の加減算：計算結果に対し、1 ULP（unit in the last place）を加算または減算して 包含性を保証する。すなわち下端からは 1 ULP を減じ、上端には 1 ULP を加える。 これは過大評価を生むが、正しさは保証される。
  • MPFR との統合：MPFR ライブラリは任意精度の浮動小数点演算において 正しい丸めを保証し、各演算で丸めモードを指定できる。 区間演算ライブラリ mpfi は MPFR をバックエンドとして利用し、 任意精度の区間演算を実現している。

区間演算の可視化

区間演算による値の包含 加法: [1.5, 2.5] + [3.0, 4.0] = [4.5, 6.5] x y x+y 1.5 2.5 3.0 4.0 + 4.5 6.5 x+y は必ずこの区間内 包含性の保証 x ∈ [a,b], y ∈ [c,d] ⇒ x+y ∈ [a+c, b+d] 依存性問題: x - x (x ∈ [0, 1]) 真の値域 {0} 区間結果 -1 1 過大評価 (overestimation)

区間演算ライブラリ

区間演算を正しく実装するには、向き付き丸めの制御や、 多数の初等関数に対する厳密な区間拡張の実装など、 相当の技術的労力を要する。 以下に主要な区間演算ライブラリを紹介する。

Arb

Fredrik Johansson により開発された球演算（ball arithmetic）ライブラリ。 区間を「中点と半径」の組 $x = (m, r)$ で表現し、$m$ は任意精度浮動小数点数、 $r$ は固定精度の正の浮動小数点数である。 すなわち、真の値は $[m - r, m + r]$ に含まれる。

この表現は対称性が良く、多くの場合に端点表現より効率的である。 Arb は複素数の球演算、多項式、行列、特殊関数など広範な機能を提供し、 SageMath や Nemo（Julia の代数システム）に統合されている。 C 言語で実装されており、GMP/MPFR/FLINT をバックエンドとして用いる。

INTLAB

Siegfried M. Rump により開発された、MATLAB/Octave 上の検証付き数値計算（verified computing）ツールボックス。 区間演算に加え、自動微分、線形・非線形方程式の検証的求解、 固有値の包含など、精度保証付き数値計算に必要な機能を網羅的に提供する。 教育・研究の両面で広く利用されている。

Boost.Interval (C++)

Boost C++ ライブラリの一部として提供される区間演算ライブラリ。 テンプレートベースの設計により、基底型（float, double など）を柔軟に選択できる。 丸めポリシーもテンプレート引数で指定可能であり、C++ での区間演算の標準的な選択肢の一つである。

filib++

高速な区間演算を目指して設計された C++ ライブラリ。 IEEE 754 の丸めモード切り替えを効率的に活用し、 基本四則演算に加えて初等関数の区間評価を高速に行う。 特に性能を重視するアプリケーションで利用される。

mpfi

MPFR をバックエンドとする多倍長精度区間演算ライブラリ。 MPFR の正しい丸めと向き付き丸めモードをそのまま活用し、 任意精度の区間演算を行うことができる。 区間の各端点は MPFR の多倍長浮動小数点数で表現される。 高精度が必要な計算（定数の厳密な包含、級数の剰余項の評価など）に適している。

ライブラリの選択指針

  • 多倍長精度が必要：Arb または mpfi
  • MATLAB/Octave 環境：INTLAB
  • C++ で高速な倍精度演算：Boost.Interval または filib++
  • 広範な数学的機能（特殊関数、多項式演算など）：Arb
  • 検証付き数値計算（線形方程式、固有値問題など）：INTLAB

まとめ

  • 歴史と動機：1960年代に Moore が提唱。浮動小数点演算の丸め誤差を自動的・厳密に追跡する枠組みである
  • 区間演算の定義：四則演算を区間に拡張し、包含性定理により真の値が結果区間に含まれることが保証される
  • 依存性問題：同一変数の複数回出現が過大評価を招く。式の書き換え、アフィン演算、Taylor模型、平均値形式が対策となる
  • 区間ニュートン法：根の存在を厳密に証明でき、不在も証明できる。Brouwer の不動点定理に基づく。二分法との組み合わせで大域的根探索が可能
  • 初等関数の区間拡張：単調関数は端点評価で済むが、冪乗 $x^n$ や三角関数は区間内の極値検出が必要。$x^2$ を専用関数として実装することで依存性による過大評価を回避できる
  • 区間行列演算：区間連立一次方程式の解集合は一般に非凸であり、Oettli-Prager 定理で特徴づけられる。最小包含区間の計算はNP困難だが、Krawczyk法やRump法で実用的に包含できる
  • 制約伝播：等式・不等式制約から区間を逆方向にも縮小する手法。大域的最適化や制約充足問題の加速に不可欠である
  • 向き付き丸め：IEEE 754 の丸めモードを活用し、浮動小数点環境で包含性を保証する。MPFR との統合が理想的
  • 主要ライブラリ：Arb（球演算）、INTLAB（MATLAB）、Boost.Interval（C++）、filib++（高速）、mpfi（多倍長精度）