線形制約を含む半正定値計画問題の解法

nAG数値計算ライブラリ > 最適化アルゴリズムExample集 > 線形制約を含む半正定値計画問題の解法

半正定値計画問題-SDP（線形制約）

このExampleでは、Pennonを用いた半正定値計画法(SDP)により、Petersen GraphのLovász theta数を求めています。Lovász theta数はグラフ理論における重要な不変量で、グラフの安定数の上界を与えます。

目的関数:

タスク	式
minimize	$x_{1}$

目的関数は $x_{1}$ の最小化で、これはLovász theta数に対応します。

$minimize x_{1}$

決定変数:

変数	範囲
$x_{i}$	$i = 1, \dots, n e + 1$

ここで、 $n e$ はグラフの辺の数を表します。よって、決定変数の数は辺の数+1となります。

制約条件:

制約	式
LMI制約	$\sum_{1 \leq i \leq j \leq n v} E_{i j} + \sum_{i = 1}^{n v} x_{1} E_{i i} + \sum_{k = 1}^{n e} x_{k + 1} (E_{v_{a} (k), v_{b} (k)} + E_{v_{b} (k), v_{a} (k)}) ⪰ O$

ここで、 $n v$ はグラフの頂点数、 $E_{i j}$ は $i$ 行 $j$ 列が1でそれ以外が0の $n v \times n v$ 行列、 $O$ は $n v \times n v$ の零行列を表します。また、 $v_{a} (k)$ と $v_{b} (k)$ は $k$ 番目の辺が接続する2つの頂点を表します。

この線形行列不等式(LMI)制約は、Lovász theta数の定義そのものになっています。

以上が、Petersen GraphのLovász theta数を半正定値計画問題として定式化した際の、決定変数、目的関数、制約条件になります。数値例としては、Petersen Graphは頂点数 $n v = 10$ 、辺数 $n e = 15$ であり、それに基づいて行列の次元等が決まります。最終的にこの問題をPennonソルバーで解くことで、Petersen GraphのLovász theta数として $4.00$ が得られます。

Exampleの実行コマンド:

python -m naginterfaces.library.examples.opt.handle_solve_pennon_lmi_ex

ソースコード表示コマンド:

python -c "import inspect; from naginterfaces.library.examples.opt import handle_solve_pennon_lmi_ex; print(''.join(inspect.getsourcelines(handle_solve_pennon_lmi_ex)[0]))"

出力結果例:

naginterfaces.library.opt.handle_solve_pennon Python Example Results.
Find the Lovasz theta number for a Petersen Graph.
Lovasz theta number of the given graph is    4.00.

マニュアル:

handle_solve_pennon_lmiのマニュアル

ソース:

#!/usr/bin/env python3
"""
``naginterfaces.library.opt.handle_solve_pennon`` Python Example,
with linear matrix inequality constraints.
"""

# nAG Copyright 2018-2019.

# pylint: disable=invalid-name,too-many-arguments,too-many-locals,too-many-statements

import numpy as np

from naginterfaces.library import opt

def main():
    """
    Example for :func:`naginterfaces.library.opt.handle_solve_pennon`,
    with linear matrix inequality constraints.

    Semidefinite programming using Pennon.

    >>> main()
    naginterfaces.library.opt.handle_solve_pennon Python Example Results.
    Find the Lovasz theta number for a Petersen Graph.
    Lovasz theta number of the given graph is    4.00.
    """

    print(
        'naginterfaces.library.opt.handle_solve_pennon Python Example Results.'
    )
    print('Find the Lovasz theta number for a Petersen Graph.')

    # The graph structure:
    # Number of vertices in the graph:
    nv = 10
    # Number of edges in the graph:
    ne = 15
    va = np.empty(ne, dtype=int)
    vb = np.empty(ne, dtype=int)
    maxe = ne
    ne = 0
    ij = [
        (1, 2),
        (2, 3),
        (3, 4),
        (4, 5),
        (1, 5),
        (1, 6),
        (2, 7),
        (3, 8),
        (4, 9),
        (5, 10),
        (6, 8),
        (6, 9),
        (7, 9),
        (7, 10),
        (8, 10),
    ]
    for idx in range(maxe):
        i = ij[idx][0]
        j = ij[idx][1]
        if 1 <= i < j <= nv:
            va[ne] = i
            vb[ne] = j
            ne += 1

    # Initial estimate of the solution:
    nvar = ne + 1
    x = [0.]*nvar

    # Create a handle for the problem:
    handle = opt.handle_init(nvar)

    # Define the quadratic objective:
    opt.handle_set_quadobj(
        handle, idxc=[1], c=[1.],
    )

    # Define the linear matrix constraint:
    nnzasum = ne + nv + (nv + 1) * nv // 2
    nnza = np.empty(nvar + 1, dtype=int)
    irowa = np.empty(nnzasum, dtype=int)
    icola = np.empty(nnzasum, dtype=int)
    a = np.empty(nnzasum)

    idx = 0
    nnza[0] = (nv + 1) * nv // 2
    for i in range(1, nv+1):
        for j in range(i, nv+1):
            irowa[idx] = i
            icola[idx] = j
            a[idx] = 1.0
            idx += 1

    nnza[1] = nv
    for i in range(1, nv+1):
        irowa[idx] = i
        icola[idx] = i
        a[idx] = 1.0
        idx += 1

    for i in range(0, ne):
        nnza[2 + i] = 1
        irowa[idx] = va[i]
        icola[idx] = vb[i]
        a[idx] = 1.0
        idx += 1

    opt.handle_set_linmatineq(
        handle, nv, nnza, irowa, icola, a,
    )

    opt.handle_opt_set(
        handle, 'Initial X = Automatic',
    )
    opt.handle_opt_set(
        handle, 'Print Level = 0',
    )
    opt.handle_opt_set(
        handle, 'Print Options = No',
    )

    # Solve the problem:
    theta = opt.handle_solve_pennon(handle, x, inform=0).rinfo[0]

    print(
        'Lovasz theta number of the given graph is {:7.2f}.'.format(theta)
    )

    # Destroy the handle:
    opt.handle_free(handle)

if __name__ == '__main__':
    import doctest
    import sys
    sys.exit(
        doctest.testmod(
            None, verbose=True, report=False,
            optionflags=doctest.ELLIPSIS | doctest.REPORT_NDIFF,
        ).failed
    )