概要

整数ベクトル $B := (b_1, \ldots, b_m) \in \mathbb{Z}^{m \times n} $ が与えられたとき，これらの整数線型結合で作られる集合 $L(B) := \{ x_1 b_1 + \ldots + x_m b_m : x_1, \ldots, x_m \in \mathbb{Z} \} \subseteq \mathbb{Z}^n $ を lattice といいます．Lattice $ L $ を生成する $ B $ を $ L $ の基底といいます．一般に $ L $ に対して複数の基底が存在します．例えば $ ( (1, 0), (0,1) ) $ と $ ( (1,2), (3,5) ) $ はどちらも $\mathbb{Z}^2$ を生成します．

Lattice における基本的な問題はできるだけシンプルな基底を求めることです．特に以下の問題たちが基本的です．

最短ベクトル問題．$L(B)$ のうち，ユークリッドノルムが最も小さいベクトルを求めよ．
最短基底問題．$L(B)$ の基底 $B'$ のうち，ユークリッドノルムが（適当な意味で）最も小さなものを求めよ．

Lenstra-Lenstra-Lovasz (1982): https://doi.org/10.1007/BF01457454 の Lattice-Basis Reduction Algorithm (LLL) はこれらの問題の近似解を求める多項式時間アルゴリズムです．近似率は $2^{\Omega(m)}$ とあまりよろしくないのですが，対象となるベクトルの本数が少なければ十分有用な解を得ることができます．

LLL の応用

無理数の有理数近似

円周率 $\pi = 3.14159265...$ を近似する有理数として $22/7 = 3.1428...$ や $355/113 = 3.14159292...$ が知られています．これらの数値を求める手法として連分数展開が知られていますが，LLL を用いて求めることもできます．

十分大きな整数 $ M $ を固定してベクトル $b_1 = (1, 0, \lfloor M \pi \rfloor)$, $b_2 = (0, 1, -M)$ を考えます．これらが生成する lattice から $(a, b, s)$ が取れたとしましょう．これは $a \lfloor M \pi \rfloor - b M = s$ を意味するため，$ a \pi - b \approx 0$ となり，$\pi$ の近似値として $b / a$ が得られました．

無理数を解にもつ代数方程式

黄金比 $ \gamma = 1.61803398... $ は代数方程式 $1 + x - x^2 = 0$ を満たすことを知っていますが，この方程式を LLL を用いて求めることができます．

十分大きな整数 $ M $ を固定してベクトル $b_1 = (1, 0, 0, M)$, $b_2 = (0, 1, 0, \lfloor M \gamma \rfloor)$, $b_3 = (0, 0, 1, \lfloor M \gamma^2 \rfloor)$ を考えます．これらが生成する lattice から $(a, b, c, s)$ が取れたとしましょう．これは $ a M + b \lfloor M \gamma \rfloor + c \lfloor M \gamma^2 \rfloor = s $ を意味するため，$ \gamma $ が $ a + b x + c x^2 \approx 0 $ を満たすことがわかりました．

小さな Bezout 係数

整数 $ a_1, \ldots, a_m $ が与えられたとき，$ g := \mathrm{gcd}(a_1, \ldots, a_m) = x_1 a_1 + \cdots + x_m a_m $ を満たす $ (x_1, \ldots, x_m) $ を Bezout 係数といいます．Bezout 係数は拡張ユークリッド互除法で求まることが知られています．Bezout 係数は一意ではありません．例えば $ \mathrm{gcd}(3141, 5926, 5358) = 1 $ について $(-1881, 997, 0)$ はひとつの Bezout 係数ですが，他にも $(81, -14, -32)$ という係数の小さな Bezout 係数が存在します．一般に，ノルム最小の Bezout 係数を求める問題は NP 困難ですが，近似でよければ LLL を用いて求められます．

記号的な無限大 $ M $ を導入して $ b_1 = (1, 0, \ldots, 0, a_1 M), \ldots, b_m = (0, \ldots, 0, 1, a_m M) $ を考えます．$ M $ は無限大なのでノルムの小さな基底において最終成分の値はできる限り小さくなるはずです．具体的にはどれか１つのベクトルだけ $g M $，残りは $0$ となります．$ g M $ があらわれるベクトルの成分は構成から Bezout 係数となり，LLL の性質から小さな値をもつことが期待できます．

整数計画問題

変数の数が定数個の整数計画問題も LLL を使うと多項式時間で解くことができます．一般に，「格子点を含まない多面体は細長い（太さが次元にのみ依存する定数になる）」ということが証明できます (Khinchin's Flatness Theorem)．多項式時間で多面体から格子点を見つけるもしくは細長い軸を求めることができるので (*)，後者の場合にその軸に沿って多面体をスライスして分割統治を行うことで整数計画問題の実行可能性問題を解くことができます．このとき分枝の個数は次元にのみ依存する定数になり，分枝の深さは次元に等しくなるので，全体として探索木のサイズが定数となります．(*) の手続きのなかで LLL が用いられます．

理論的に価値の大きな結果ですが，実用的には通常の分枝限定法のほうが遥かに早いと言われています．

LLL の実装

LLL は基底に対して適当な基準で行基本変形を繰り返すタイプのアルゴリズムです．基準をうまく選ぶことによって手続き中でノルムが減少していくことが示せます．正当性の証明は複雑ですが，実装自体は Gauss の消去法やら拡張互除法やらとそれほど大きく違わないシンプルなものとなります．

Wikipedia.en: https://en.wikipedia.org/wiki/LLL_algorithm に擬似コードがあります．これは Lenstra-Lenstra-Lovasz によるオリジナル実装に近いものですが，有理数を管理しなければいけないため定数倍がよくありません．有理数を必要としない実装が de Weger (1987): https://doi.org/10.1016/0022-314X(87)90088-6 によって提案されています．以下はそれを実装したものです．

注意として，LLL を有効に使うには原則として多倍長演算が必要です．これは行基本変形中に係数が増大するためですが，それでも高々多項式程度に収まるのが LLL の偉いところです．

use num::traits::{NumRef, RefNum};
use num::{Integer, Signed};

pub trait FracRound {
    fn frac_round(&self, other: &Self) -> Self;
}
impl<T> FracRound for T
where
    T: Clone + Integer + Signed + std::fmt::Debug,
{
    fn frac_round(&self, other: &Self) -> Self {
        let sign = self.signum() * other.signum();
        let num = self.abs();
        let den = other.abs();
        let (div, rem) = num.div_rem(&den);
        let value = match (rem.clone() + rem.clone()).cmp(&den) {
            std::cmp::Ordering::Less => div,
            std::cmp::Ordering::Equal => {
                if sign.is_positive() {
                    div
                } else {
                    div + T::one()
                }
            }
            std::cmp::Ordering::Greater => div + T::one(),
        };
        sign * value
    }
}

pub struct LatticeBasisReductionDeWeger<T> {
    b: Vec<Vec<T>>,
    u: Vec<Vec<T>>,
    lambda: Vec<Vec<T>>,
    d: Vec<T>,
}
impl<T> LatticeBasisReductionDeWeger<T>
where
    T: Clone + Integer + Signed + NumRef + From<i8> + std::fmt::Debug,
    for<'a> &'a T: RefNum<T>,
{
    fn dot(a: &[T], b: &[T]) -> T {
        a.iter()
            .zip(b.iter())
            .fold(T::zero(), |sum, (a, b)| sum + a * b)
    }

    pub fn new<S: Clone + Into<T>>(matrix: Vec<Vec<S>>) -> Self {
        let m = matrix.len();
        let mut b = Vec::new();
        for row in matrix {
            b.push(row.into_iter().map(|x| x.into()).collect::<Vec<_>>());
        }
        let mut c = b.clone();
        let mut u = vec![vec![T::zero(); m]; m];
        for i in 0..m {
            u[i][i] = T::one();
        }
        let mut d = vec![T::one(); m];
        let mut lambda = vec![vec![T::zero(); m]; m];
        for i in 0..m {
            for j in 0..i {
                lambda[i][j] = Self::dot(&b[i], &c[j]);
                for p in 0..c[i].len() {
                    c[i][p] = &d[j] * &c[i][p] - &lambda[i][j] * &c[j][p];
                    if j >= 1 {
                        c[i][p] = &c[i][p] / &d[j - 1];
                    }
                }
            }
            d[i] = Self::dot(&c[i], &c[i]);
            if i >= 1 {
                d[i] = &d[i] / &d[i - 1];
            }
        }

        Self { b, u, lambda, d }
    }

    fn reduce(&mut self, k: usize, l: usize) {
        if T::from(2) * self.lambda[k][l].abs() > self.d[l] {
            let r = self.lambda[k][l].frac_round(&self.d[l]);

            for p in 0..self.b[k].len() {
                self.b[k][p] = &self.b[k][p] - &r * &self.b[l][p];
            }
            for p in 0..self.u[k].len() {
                self.u[k][p] = &self.u[k][p] - &r * &self.u[l][p];
            }
            for j in 0..l {
                self.lambda[k][j] = &self.lambda[k][j] - &r * &self.lambda[l][j];
            }
            self.lambda[k][l] = &self.lambda[k][l] - &r * &self.d[l];
        }
    }

    fn swap(&mut self, k: usize) {
        self.b.swap(k, k - 1);
        self.u.swap(k, k - 1);
        for j in 0..k - 1 {
            // swap lambda[k-1][j] and lambda[k][j]
            let (a, b) = self.lambda.split_at_mut(k);
            std::mem::swap(&mut a[k - 1][j], &mut b[0][j]);
        }
        for i in k + 1..self.lambda.len() {
            let t =
                &self.lambda[i][k - 1] * &self.d[k] - &self.lambda[i][k] * &self.lambda[k][k - 1];
            self.lambda[i][k - 1] = (&self.lambda[i][k - 1] * &self.lambda[k][k - 1]
                + if k >= 2 {
                    &self.lambda[i][k] * &self.d[k - 2]
                } else {
                    self.lambda[i][k].clone()
                })
                / &self.d[k - 1];
            self.lambda[i][k] = t / &self.d[k - 1];
        }
        self.d[k - 1] = (if k >= 2 {
            &self.d[k - 2] * &self.d[k]
        } else {
            self.d[k].clone()
        } + &self.lambda[k][k - 1] * &self.lambda[k][k - 1])
            / &self.d[k - 1];
    }

    pub fn run(&mut self) {
        let m = self.b.len();
        let mut k = 1;
        while k < m {
            self.reduce(k, k - 1);

            let lhs = if k >= 2 {
                &self.d[k - 2] * &self.d[k] - &self.lambda[k][k - 1] * &self.lambda[k][k - 1]
            } else {
                &self.d[k] - &self.lambda[k][k - 1] * &self.lambda[k][k - 1]
            };
            let rhs = &self.d[k - 1] * &self.d[k - 1];

            if T::from(4) * lhs < T::from(3) * rhs {
                self.swap(k);
                if k >= 2 {
                    k -= 1;
                }
            } else {
                for i in (0..k - 1).rev() {
                    self.reduce(k, i);
                }
                k += 1;
            }
        }
    }
}

TM's memo

Lenstra-Lenstra-Lovasz's Lattice Basis Reduction

概要