Processing math: 100%

半正定値計画法(SDP)を用いた最近接相関行列

nAG Library for Python Example集

このページは、nAGライブラリのJupyterノートブックExampleの日本語翻訳版です。オリジナルのノートブックはインタラクティブに操作することができます。

from naginterfaces.library import opt, lapackeig, correg
from naginterfaces.base import utils
import numpy as np

半正定値計画法(SDP)を用いた最近接相関行列

このノートブックの正しいレンダリング

このノートブックでは、方程式と参照のためにlatex_envs Jupyter拡張機能を使用しています。LaTeXがローカルのJupyterインストールで適切にレンダリングされない場合、この拡張機能をインストールしていない可能性があります。詳細は https://jupyter-contrib-nbextensions.readthedocs.io/en/latest/nbextensions/latex_envs/README.html を参照してください。

nAGライブラリのインストールとこのノートブックの実行

このノートブックを実行するにはnAGライブラリ for Pythonが必要です。ライブラリのダウンロード、インストール、ライセンス取得については、index.htmファイルの指示をお読みください。

ノートブックの実行方法についてはこちらをご覧ください。

はじめに

相関行列ではない行列Gから始めます: [1100111001110012]

確かに対称で対角線上に1がありますが、固有値を見てみると:

G = np.array([[1, -1,  0, 0],
              [-1, 1, -1, 0],
              [0, -1, 1, -1],
              [0, 0, -1, 1]],
             dtype=float)
n = G.shape[0]

# Gの固有値を計算する
itype = 1
jobz = 'N'
uplo = 'U'
A = G.copy()
B = np.identity(n, dtype=float)
_, _, sv = lapackeig.dsygv(itype, jobz, uplo, A, B)
print('Minimum eigenvalue:', min(sv))
Minimum eigenvalue: -0.6180339887498952

負の固有値があります!したがって、Gは半正定値ではなく、相関行列ではありません。

Gの摂動

Gに摂動を加えて、フロベニウスノルムの意味でGに最も近い相関行列G*を得ることができます: min||GG||Fs.t G0 ここで、 G=[11+x1001+x111+x2001+x211+x3001+x32] これはSDPの定義に完全に適合します。まず、変数の数を計算して定義しましょう:

# 変数の数: 対角線より上の非ゼロ要素の数
nvar = 0
for i in range(0, n):
    for j in range(i+1, n):
        if G[i, j] != 0.0:
            nvar += 1

# 問題ハンドルを初期化する
handle = opt.handle_init(nvar)

次に目的関数を設定します。GGのフロベニウスノルムは単純に以下のようになります: ix2i

# 私たちの変数はベクトルとして格納されているため、単に最小化します
# 修正の二乗和 --> H は単位行列、c = 0。
h = np.ones(nvar, dtype=float)
irowh = np.arange(nvar, dtype=int) + 1
icolh = np.arange(nvar, dtype=int) + 1
opt.handle_set_quadobj(handle, idxc=None, c=None, irowh=irowh, icolh=icolh, h=h)

ここで制約 G0 を定義する必要があります。これは線形行列不等式で行うことができます: x1[11]+x2[11]+x3[11]+G0 対称行列のみを扱うため、上三角部分のみを定義する必要があります

# 正定値制約: 線形行列不等式
# 非ゼロ要素の数の合計 => Gの完全三角形
# 1 各E_ijについて
nnzasum = int(n*(n+1)/2 + nvar)

# 線形結合 G + x_ijE_ij における各行列の非ゼロ要素の数
nnza = np.empty(nvar+1, dtype=int)
nnza[0] = n*(n+1)/2
nnza[1:nvar+1] = 1

# 行列不等式を3つの疎配列に格納する
irowa = np.empty(nnzasum, dtype=int)
icola = np.empty(nnzasum, dtype=int)
a =  np.empty(nnzasum, dtype=float)
# G の上三角部分を -A_0 にコピーする
idx = 0
for i in range(n):
    for j in range(i, n):
        irowa[idx] = i + 1
        icola[idx] = j + 1
        a[idx] = -G[i, j]
        idx += 1
# E_ijは1つの非ゼロ要素を持つ
for i in range(n):
    for j in range(i+1, n):
        if (G[i,j] != 0):
            irowa[idx] = i + 1
            icola[idx] = j + 1
            a[idx] = 1.0
            idx += 1

dima = n
idblk = opt.handle_set_linmatineq(handle, dima, nnza, irowa, icola, a, blksizea=None, idblk=0)

問題が完全に定義されたので、ソルバーに渡すことができます

# 入出力
iom = utils.FileObjManager(locus_in_output=False)

# オプション引数を設定する
for option in ['Print Options = No',
               'Initial X = Automatic']:
    opt.handle_opt_set(handle, option)
    
x = np.empty(nvar)
inform = 0
x, _, _, _, rinfo, stats, _ = opt.handle_solve_pennon(handle, x, inform, u=None, uc=None, ua=None, io_manager=iom)
 --------------------------------
  E04SV, NLP-SDP Solver (Pennon)
 --------------------------------

 Problem Statistics
   No of variables                  3
     bounds               not defined
   No of lin. constraints           0
     nonzeroes                      0
   No of matrix inequal.            1
     detected matrix inq.           1
       linear                       1
       nonlinear                    0
       max. dimension               4
     detected normal inq.           0
       linear                       0
       nonlinear                    0
   Objective function       Quadratic

 --------------------------------------------------------------
  it|  objective |  optim  |   feas  |  compl  | pen min |inner
 --------------------------------------------------------------
   0  0.00000E+00  0.00E+00  6.19E-01  6.63E+00  1.00E+00   0
   1  4.12017E-01  6.38E-04  0.00E+00  1.44E+00  1.00E+00   5
   2  3.29642E-01  7.76E-04  0.00E+00  4.96E-01  4.65E-01   2
   3  2.65315E-01  1.02E-04  0.00E+00  1.55E-01  2.16E-01   3
   4  2.33229E-01  1.03E-03  0.00E+00  4.71E-02  1.01E-01   3
   5  2.19082E-01  2.22E-03  0.00E+00  1.46E-02  4.68E-02   3
   6  2.13121E-01  2.12E-03  0.00E+00  4.72E-03  2.18E-02   3
   7  2.10698E-01  1.26E-03  0.00E+00  1.56E-03  1.01E-02   3
   8  2.09756E-01  4.90E-04  0.00E+00  4.85E-04  4.71E-03   3
   9  2.09413E-01  1.13E-04  0.00E+00  1.21E-04  2.19E-03   3
  10  2.09310E-01  1.95E-03  0.00E+00  1.63E-05  1.02E-03   2
  11  2.09297E-01  1.25E-05  0.00E+00  2.77E-06  4.74E-04   2
  12  2.09294E-01  2.68E-07  0.00E+00  3.89E-07  2.21E-04   2
  13  2.09294E-01  2.25E-09  0.00E+00  5.43E-08  1.03E-04   2
 --------------------------------------------------------------
 Status: converged, an optimal solution found
 --------------------------------------------------------------
 Final objective value                2.092940E-01
 Relative precision                   2.759238E-07
 Optimality                           2.249025E-09
 Feasibility                          0.000000E+00
 Complementarity                      5.426796E-08
 Iteration counts
   Outer iterations                             13
   Inner iterations                             36
   Linesearch steps                             36
 Evaluation counts
   Augm. Lagr. values                           50
   Augm. Lagr. gradient                         50
   Augm. Lagr. hessian                          36
 --------------------------------------------------------------
# 最も近い相関行列を形成する
G_star = np.zeros((n,n))
idx = 0
for i in range(n):
    for j in range (i+1, n):
        if G[i,j] != 0:
            G_star[i,j] = G[i,j] + x[idx]
            G_star[j,i] += G[i,j] + x[idx]
            idx += 1
    G_star[i,i] = 1.0
print('The Nearest correlation matrix to G is')
print(G_star)

Gに最も近い相関行列は [[ 1. -0.68232776 0. 0. ] [-0.68232776 1. -0.53442871 0. ] [ 0. -0.53442871 1. -0.68232776] [ 0. 0. -0.68232776 1. ]]

ここで、その固有値をテストしてみましょう

# Gの固有値を計算する
itype = 1
jobz = 'N'
uplo = 'U'
A = G_star.copy()
B = np.identity(n, dtype=float)
_, _, sv = lapackeig.dsygv(itype, jobz, uplo, A, B)
print('Minimum eigenvalue:', min(sv))
Minimum eigenvalue: 7.946753527371534e-08

専用のソルバーが存在します!

この問題は十分に一般的であるため、nAGライブラリで専用のソルバーが開発されました。試してみましょう!

G_star, _, _, _ = correg.corrmat_nearest(G)
print(G_star)

[[ 1. -0.8084125 0.1915875 0.10677505] [-0.8084125 1. -0.65623269 0.1915875 ] [ 0.1915875 -0.65623269 1. -0.8084125 ] [ 0.10677505 0.1915875 -0.8084125 1. ]]

しかし、疎性パターンは保持されていません。別のソルバーを使用すると、一部の値を固定することができます:

# G_star[i, j] は H[i,j]==1 の場合、G[i,j] に固定されます
H = np.array([[0, 0, 1, 1],
             [0, 0, 0, 1],
             [1, 0, 0, 0],
             [1, 1, 0, 0]])
alpha = 0.0
G_star, _, _ = correg.corrmat_fixed(G, alpha, H, 4)
print(G_star)

[[ 1. -0.6823278 0. 0. ] [-0.6823278 1. -0.53442877 0. ] [ 0. -0.53442877 1. -0.6823278 ] [ 0. 0. -0.6823278 1. ]]

help(opt.handle_solve_pennon)

以下は翻訳結果です:

naginterfaces.libraryモジュールのhandle_solve_pennon関数に関するヘルプ:

handle_solve_pennon(handle, x, inform, u=None, uc=None, ua=None, io_manager=None)
    ``handle_init``で初期化され、半正定値計画法(SDP)や双線形行列不等式(BMI)を含む
    SDPなど、スイートの他の関数で定義された互換性のある問題に対して
    Pennonソルバーを実行します。
    
    注意: この関数はオプションのアルゴリズムパラメータを使用します。以下も参照してください:
    ``handle_opt_set``、``handle_opt_get``。
    
    ``handle_solve_pennon``は、二次計画法(QP)、線形半正定値計画法(SDP)、
    双線形行列不等式を含む半正定値計画法(BMI-SDP)などの問題に対する
    nAG最適化モデリングスイートのソルバーです。
    
    詳細な情報については、e04svに関するnAGライブラリドキュメントを参照してください
    
    https://www.nag.com/numeric/nl/nagdoc_27.1/flhtml/e04/e04svf.html
    
    パラメータ
    ----------
    handle : Handle
        問題へのハンドル。初期化されている必要があり(例えば``handle_init``によって)、
        ``handle_solve_pennon``と互換性のある問題定式化を保持している必要があります。
        nAG最適化モデリングスイートの呼び出し間で変更してはいけません。
    
    x : float, array-like, shape (nvar)
        注意: 中間停止は'Monitor Frequency' > 0の場合にのみ発生します。
    
        'Initial X' = 'USER'(デフォルト)の場合、x^0は変数xの初期推定値です。
        それ以外の場合、`x`を設定する必要はありません。
    
        `中間エントリー時`: 入力は無視されます。
    
    inform : int
        注意: 中間停止は'Monitor Frequency' > 0の場合にのみ発生します。
    
        `初期エントリー時`: 効果はありません。
    
        `中間エントリー時`: 0に設定されると、現在の問題の解決が終了し、
        関数は`errno` = 20を返します。それ以外の場合、関数は続行します。
    
    u : None または float, array-like, shape (nnzu), オプション
        注意: 中間停止は'Monitor Frequency' > 0の場合にのみ発生します。
    
        注意: nnzu > 0の場合、`u`は(標準的な)制約のラグランジュ乗数(双対変数)を
        保持します。つまり、``handle_set_simplebounds``で定義された単純な境界と、
        ``handle_set_linconstr``で定義されたm_B個の線形制約のセットです。
        初期推定値、中間近似値、または最終値のいずれかです。
        [ラグランジュ乗数の構造]を参照してください。
    
        注意: nnzu = 0の場合、`u`は参照されません。
    
        'Initial U' = 'USER'の場合(デフォルトは'AUTOMATIC')、u^0はラグランジュ乗数u
        の初期推定値です。それ以外の場合、`u`を設定する必要はありません。
    
        `中間エントリー時`: 入力は無視されます。
    
    uc : None または float, array-like, shape (nnzuc), オプション
        `uc`は、二次錐制約の定義を可能にするnAGライブラリの将来のリリースのために
        予約されています。
    
        現時点では参照されません。
    
    ua : None または float, array-like, shape (nnzua), オプション
        注意: 中間停止は'Monitor Frequency' > 0の場合にのみ発生します。
    
        nnzua > 0の場合、`ua`は``handle_set_linmatineq``で定義された行列制約の
        ラグランジュ乗数を保持します。
       ``handle_set_quadmatineq``については、[ラグランジュ乗数の構造]を参照してください。

        nnzua = 0の場合、`ua`は参照されません。

        'Initial U' = 'USER'の場合(デフォルトは'AUTOMATIC')、U^0は行列ラグランジュ乗数Uの初期推定値です。
        それ以外の場合、`ua`を設定する必要はありません。

        `中間エントリー時`: 入力は無視されます。

    io_manager : FileObjManager, オプション
        このルーチンでのI/O用マネージャー。

    戻り値
    -------
    x : float, ndarray, shape (nvar)
        `中間終了時`: 現在の外部反復の終了時における変数xの値。

        `最終終了時`: 変数xの最終値。

    u : float, ndarray, shape (nnzu)
        `中間終了時`: 現在の外部反復の終了時における乗数uの推定値。乗数uの最終値。

    uc : float, ndarray, shape (nnzuc)
        `uc`は、二次錐制約の定義を可能にするnAGライブラリの将来のリリースのために予約されています。

        現時点では参照されません。

    ua : float, ndarray, shape (nnzua)
        `中間終了時`: 外部反復の終了時における行列乗数Uの推定値。

        `最終終了時`: 乗数Uの最終推定値。

    rinfo : float, ndarray, shape (32)
        `中間または最終エントリー時`: 現在の(または最終の)外部反復の終了時における誤差測定値と様々な指標(詳細は[アルゴリズムの詳細]を参照)。以下の表に示されています:

        [表省略]

    stats : float, ndarray, shape (32)
        `中間または最終終了時`: 現在の(または最終の)外部反復の終了時におけるソルバー統計。以下の表に示されています。時間統計は'Stats Time'が設定されている場合のみ提供され(デフォルトは'NO')、測定時間は秒単位で返されます。

        [表省略]

    inform : int
        注:中間停止は'Monitor Frequency' > 0の場合にのみ発生します。

        `中間終了時`: `inform` = 1。

        `最終終了時`: `inform` = 0。

    その他のパラメータ
    ----------------
    'Defaults' : 値なし
        この特別なキーワードを使用して、すべてのオプションをデフォルト値にリセットできます。
        このキーワードで与えられた値は無視されます。

    'DIMACS Measures' : str
        デフォルト = 'CHECK'

        問題が線形半正定値計画問題の場合、この引数はDIMACS誤差測定([停止基準]参照)を計算および/またはチェックするかどうかを指定します。
        その他の場合、このオプションは自動的に'NO'に戻ります。

        制約:'DIMACS Measures' = 'COMPUTE'、'CHECK'または'NO'。

    'Hessian Density' : str
        デフォルト = 'AUTO'

        このオプションは、拡張ラグランジアンF(x,u,v,U,p,P)のヘッセ行列をどのように構築するかについてソルバーに指示します。
        'AUTO'オプションは決定をソルバーに任せ、推奨オプションです。
        'DENSE'に設定すると、自動検出をバイパスし、ヘッセ行列は
       は常に密行列として構築されます。
        オプション 'SPARSE' は、可能な場合、ソルバーに行列のスパース格納と
        因数分解を使用するよう指示します。
    
        制約: 'Hessian Density' = 'AUTO'、'DENSE' または 'SPARSE'
    
    'Infinite Bound Size' : float
        デフォルト = 10^20
    
        これは問題制約の定義における「無限」境界 bigbnd を定義します。
        bigbnd 以上の上界は +無限大とみなされます(同様に、
        -bigbnd 以下の下界は -無限大とみなされます)。
        このオプションの変更は、既に定義された制約には影響しません。
        変更後に定式化された制約のみが影響を受けます。
    
        制約: 'Infinite Bound Size' >= 1000.
    
    'Initial P' : str
        デフォルト = 'AUTOMATIC'
    
        このオプションは、ペナルティオプション p^0、P^0 の選択を定義します。
        [アルゴリズム1]を参照してください。
    
        'Initial P' = 'AUTOMATIC'
            ペナルティオプションは、オプション 'Init Value P'、'Init Value Pmat' に
            設定され、自動スケーリングの対象として自動的に選択されます。
            開始点でのすべての行列制約に対してペナルティ関数 Phi_P() が
            定義されるように、P^0 が増加する可能性があることに注意してください。
    
        'Initial P' = 'KEEP PREVIOUS'
            可能であれば、ペナルティオプションはソルバーの前回の実行から
            保持されます。そうでない場合、このオプションは 'AUTOMATIC' に
            戻ります。行列ペナルティオプションが前回の実行と同じであっても、
            開始点でペナルティ関数 Phi_P() が適切に定義されるように
            増加の対象となる可能性があることに注意してください。
    
        制約: 'Initial P' = 'AUTOMATIC' または 'KEEP PREVIOUS'
    
    'Initial U' : str
        デフォルト = 'AUTOMATIC'
    
        この引数は、どの初期ラグランジュ乗数を使用するかについて
        ソルバーに指示します。
    
        'Initial U' = 'AUTOMATIC'
            ラグランジュ乗数は自動スケーリングによって自動的に選択されます。
    
        'Initial U' = 'USER'
            配列 `u` と `ua`(提供されている場合)の値が、自動調整の対象となる
            初期ラグランジュ乗数として使用されます。いずれかの配列が提供
            されていない場合、欠落データの選択は 'AUTOMATIC' と同様です。
    
        'Initial U' = 'KEEP PREVIOUS'
            ラグランジュ乗数はソルバーの前回の実行から保持されます。
            このオプションが最初の実行に設定されているか、オプションが
            ソルバーのアプローチを変更する場合、選択は自動的に 'AUTOMATIC' に
            戻ります。これは、例えば、解の精度を高めるためにソルバーを
            ホットスタートする場合に有用かもしれません。
    
        制約: 'Initial U' = 'AUTOMATIC'、'USER' または 'KEEP PREVIOUS'
    
    'Initial X' : str
        デフォルト = 'USER'
    
        この引数は、どの開始点 x^0 を使用するかを指示します。
    
        'Initial X' = 'AUTOMATIC'
            開始点は、変数の単純な境界を満たすか、またはゼロとして
            自動的に選択されます
           ベクトル。引数 `x` の入力は無視されます。

    'Initial X' = 'USER'
            引数 `x` の初期値が開始点として使用されます。

    制約: 'Initial X' = 'AUTOMATIC' または 'USER'。

    'Init Value P' : float
        デフォルト = 1.0

        この引数は、(標準の)不等式に対する初期ペナルティオプション p^0 の値を定義します。
        ペナルティの低い値は、内部問題の解を実行可能領域に、つまり望ましい結果により近づけます。
        しかし、それはまたシステムの不良条件を増加させます。
        良い開始点が提供されない限り、ペナルティを低く設定することはお勧めできません。

        制約: fourthroot(epsilon) <= 'Init Value P' <= 10^4。

    'Init Value Pmat' : float
        デフォルト = 1.0

        このオプションの値は、行列不等式に対する初期ペナルティオプション P^0 を提案します。
        これは 'Init Value P' に似ています(そして同じアドバイスが適用されます)が、
        行列制約が実際のペナルティオプションよりも実行不可能な場合、P^0 は自動的に増加します。

        制約: fourthroot(epsilon) <= 'Init Value Pmat' <= 10^4。

    'Inner Iteration Limit' : int
        デフォルト = 100

        各外部反復で [アルゴリズム2] が実行する内部反復(ニュートンステップ)の最大数。
        このオプションを低く設定しすぎると、`errno` = 23 につながる可能性があります。
        100を超える値は収束を改善する可能性が低いです。

        制約: 'Inner Iteration Limit' > 0。

    'Inner Stop Criteria' : str
        デフォルト = 'HEURISTIC'

        内部サブ問題の解の精度 alpha は [アルゴリズム1] で決定され、通常の状況では
        [アルゴリズム2] が与えられた 'Inner Iteration Limit' 内でこの精度に達することが期待されます。
        問題が検出され、'Inner Stop Criteria' = 'HEURISTIC' の場合、[アルゴリズム2] は
        要求された精度または 'Inner Iteration Limit' に達する前に停止することが許可されます。
        これは通常、多くの無駄な反復を節約し、ソルバーは次の反復で回復する可能性があります。
        ヒューリスティックな問題検出があなたの問題に適していないと疑う場合、
        'Inner Stop Criteria' = 'STRICT' を設定することでそのような動作を禁止します。

        制約: 'Inner Stop Criteria' = 'HEURISTIC' または 'STRICT'。

    'Inner Stop Tolerance' : float
        デフォルト = 10^-2

        このオプションは、[アルゴリズム2] によって解かれる最初の内部問題に対して
        必要な精度 alpha^0 を設定します。
        内部問題の解の精度は、最初の外部反復では非常に高である必要はなく、
        最後の反復で最適性の限界 epsilon_2 に達するように外部反復を通じて自動的に調整されます。

        alpha^0 を制限的すぎる(低すぎる)に設定すると、最初の外部反復で必要な
        内部反復の数が増加し、`errno` = 23 につながる可能性があります。
        特定のケースでは、より緩和された(高い)alpha^0 を使用し、'P Update Speed' を
        増加させることが有用な場合があります。これにより、開始時に必要な内部反復の数が減少するはずです。
       外部反復の回数が増える可能性と引き換えに計算を簡略化します。

    制約: epsilon < '内部停止許容値' <= 10^3。

    'ラインサーチモード' : str
        デフォルト = 'AUTO'

        これは[アルゴリズム2]のステップサイズ選択を制御します。
        'ラインサーチモード' = 'FULLSTEP'(線形問題のデフォルト)の場合、
        可能な限り単位ステップが取られ、ステップの短縮は
        行列ペナルティ関数Phi_P()((6)参照)の未定義領域を
        避けるためにのみ行われます。
        これは線形問題には使用できますが、非線形問題には
        推奨されません。
        'ラインサーチモード' = 'ARMIJO'の場合、代わりにArmijoの
        バックトラッキングラインサーチが使用されます。これはかなり
        基本的なラインサーチです。
        'ラインサーチモード' = 'GOLDSTEIN'の場合、Goldstein条件に
        基づく3次のセーフガードされたラインサーチが採用されます。
        これが非線形問題に推奨される(そしてデフォルトの)選択肢です。

        制約: 'ラインサーチモード' = 'AUTO'、'FULLSTEP'、'ARMIJO'または
        'GOLDSTEIN'。

    'リスト' : str
        デフォルト = 'NO'

        各オプション指定の印刷をオンにしたい場合、
        この引数を'YES'に設定することができます。

        制約: 'リスト' = 'YES'または'NO'

    'モニター頻度' : int
        デフォルト = 0

        'モニター頻度' > 0の場合、ソルバーはi回目の
        外部反復の終わりごとに戻ってきます。
        これらの中間終了時には、現在の点`x`と
        ラグランジュ乗数`u`、`ua`(要求された場合)が
        統計情報とエラー測定値(`rinfo`、`stats`)とともに
        提供されます。
        引数`inform`は中間終了と最終終了を区別し、
        即時終了も可能にします。

        'モニター頻度' = 0の場合、ソルバーは最終点で
        一度だけ停止し、中間終了は行われません。

        制約: 'モニター頻度' >= 0。

    'モニタリングファイル' : int
        デフォルト = -1

        i >= 0の場合、二次(モニタリング)出力のユニット番号。
        -1に設定すると、二次出力は提供されません。
        以下の情報がユニットに出力されます:

        -   オプションのリスト;

        -   問題の統計、反復ログ、および'モニタリングレベル'で
            設定された最終状態。

        制約: 'モニタリングファイル' >= -1。

    'モニタリングレベル' : int
        デフォルト = 4

        この引数は、ソルバーが二次出力に印刷する
        情報の詳細量を設定します。
        レベルの意味は'印刷レベル'と同じです。

        制約: 0 <= 'モニタリングレベル' <= 5。

    '外部反復制限' : int
        デフォルト = 100

        [アルゴリズム1]が実行する外部反復の最大数。
        '外部反復制限' = 0の場合、反復は実行されず、
        停止基準に必要な量のみが計算され、`rinfo`で返されます。
        これは'初期X' = 'ユーザー'および'初期U' = 'ユーザー'と
        関連して、与えられた点の最適性をチェックするのに
        有用かもしれません。
        ただし、与えられた点の可能な修正に関するルールは
       開始点はまだ適用されます。`u`と`ua`を参照してください。
        オプションを低すぎに設定すると、`errno` = 22 になる可能性があります。
    
        制約: '外部反復回数制限' >= 0
    
    'P 最小値' : float
        デフォルト = sqrt(epsilon)
    
        これは、(標準的な)不等式に使用される最小のペナルティ値pであるp_minを制御します。
        一般に、ペナルティオプションの値が非常に小さいと、不適切な条件が発生し、数値的な困難を引き起こす可能性があります。
        一方、p_minが非常に高いと、アルゴリズムが要求された実行可能性の精度に到達するのを妨げます。
        通常の状況下では、デフォルト値が推奨されます。
    
        制約: epsilon <= 'P 最小値' <= 10^-2
    
    'Pmat 最小値' : float
        デフォルト = sqrt(epsilon)
    
        これは、最小行列ペナルティオプションP_minに対する'P 最小値'の等価物です。
        同じアドバイスが適用されます。
    
        制約: epsilon <= 'Pmat 最小値' <= 10^-2
    
    '優先度' : str
        デフォルト = 'SPEED'
    
        このオプションは、行列制約(6)からシステムヘッセ行列への寄与がどのように計算されるかに影響します。
        デフォルトオプションの'優先度' = 'SPEED'は、ほとんどの場合に適しているはずです。
        ただし、非常に高次元の行列制約を扱う場合、顕著なメモリオーバーヘッドが発生する可能性があり、
        '優先度' = 'MEMORY'に切り替える必要がある場合があります。
    
        制約: '優先度' = 'SPEED' または 'MEMORY'
    
    '前処理ブロック検出' : str
        デフォルト = 'YES'
    
        '前処理ブロック検出' = 'YES'の場合、前処理中に行列制約がチェックされ、
        より小さな独立したものに分割できるかどうかが判断され、ソルバーの速度が向上します。
    
        制約: '前処理ブロック検出' = 'YES' または 'NO'
    
    '出力ファイル' : int
        デフォルト = アドバイザリメッセージユニット番号
    
        i >= 0の場合、ソルバーの主要出力のユニット番号です。
        '出力ファイル' = -1の場合、他の設定に関係なく主要出力が完全にオフになります。
        デフォルト値は、オプションの初期化時、例えばハンドルの初期化時のアドバイザリメッセージユニット番号です。
        以下の情報がユニットに出力されます:
    
        -   '出力オプション'で設定された場合のオプションのリスト;
    
        -   '出力レベル'で設定された問題の統計、反復ログ、およびソルバーからの最終ステータス。
    
        制約: '出力ファイル' >= -1
    
    '出力レベル' : int
        デフォルト = 2
    
        この引数は、ソルバーが主要出力にどの程度詳細な情報を出力するかを定義します。
    
        [表省略]
    
        制約: 0 <= '出力レベル' <= 5
    
    '出力オプション' : str
        デフォルト = 'YES'
    
        '出力オプション' = 'YES'の場合、オプションのリストが主要出力に印刷されます。
    
        制約: '出力オプション' = 'YES' または 'NO'
    
    'P 更新速度' : int
        デフォルト = 12
    
        このオプションは、ペナルティオプションp,Pが更新される速度([アルゴリズム1]、ステップ(3))に影響し、
        間接的に外部反復の全体数に影響を与えます。
       その値は、初期ペナルティ値p^0、P^0から最小値p_minとP_minの半分に到達するのに必要な典型的な外部反復回数として解釈できます。
        3未満の値は非常に積極的なペナルティ更新戦略を引き起こし、内部反復回数の増加や数値的な困難につながる可能性があります。
        一方、15を超える値は比較的保守的なアプローチを生み出し、外部反復回数の増加につながります。
    
        ソルバーが問題で困難に遭遇した場合、より高い値が役立つかもしれません。
        問題がうまく動作している場合、より低い値を設定すると速度が向上する可能性があります。
    
        制約: 1 <= 'P Update Speed' <= 100.
    
    'Stats Time' : str
        デフォルト = 'NO'
    
        この引数は、アルゴリズムの様々な部分のタイミングをオンにして、時間の大部分がどこで費やされているかをより良く把握できるようにします。
        これは異なる解決アプローチの選択に役立つかもしれません。
        CPUとウォールクロック時間のどちらかを選択することができます。
        'YES'の選択は'WALL CLOCK'と同等です。
    
        制約: 'Stats Time' = 'YES'、'NO'、'CPU'または'WALL CLOCK'。
    
    'Stop Criteria' : str
        デフォルト = 'SOFT'
    
        'Stop Criteria' = 'SOFT'の場合、ソルバーは、より良い解の推定に到達できないと予測した場合、`errno` = 50で最適でない解で早期に停止することが許可されます。
        これが推奨オプションです。
    
        制約: 'Stop Criteria' = 'SOFT'または'STRICT'。
    
    'Stop Tolerance 1' : float
        デフォルト = max(10^-6, sqrt(epsilon))
    
        このオプションは、相対双対ギャップ(0)と相対精度(1)の許容誤差として使用されるepsilon_1を定義します。[停止基準]を参照してください。
    
        制約: 'Stop Tolerance 1' > epsilon.
    
    'Stop Tolerance 2' : float
        デフォルト = max(10^-7, sqrt(epsilon))
    
        このオプションは、KKT条件からの最適性(2)と相補性(4)のテスト、または'DIMACS Measures' = 'Check'の場合はすべてのDIMACSエラー測定に代わって使用されるepsilon_2の値を設定します。
        [停止基準]を参照してください。
    
        制約: 'Stop Tolerance 2' > epsilon.
    
    'Stop Tolerance Feasibility' : float
        デフォルト = max(10^-7, sqrt(epsilon))
    
        この引数は、解の実行可能性(3)、epsilon_feasの受け入れ限界を設定します。
        [停止基準]を参照してください。
    
        制約: 'Stop Tolerance Feasibility' > epsilon.
    
    'Task' : str
        デフォルト = 'MINIMIZE'
    
        この引数は、最適化の必要な方向を指定します。
        'Task' = 'FEASIBLE POINT'の場合、目的関数(設定されている場合)は無視され、与えられた許容誤差に関して実行可能な点が見つかるとすぐにアルゴリズムは停止します。
        目的関数が設定されていない場合、'Task'は自動的に'FEASIBLE POINT'に戻ります。
    
        制約: 'Task' = 'MINIMIZE'、'MAXIMIZE'または'FEASIBLE POINT'。
    
    'Transform Constraints' : str
        デフォルト = 'AUTO'
    
        この引数は、等式制約がソルバーによってどのように扱われるかを制御します
       ソルバー。
        'Transform Constraints' = 'EQUALITIES'の場合、(4)のすべての等式
        制約h_k(x) = 0は、2つの不等式
        h_k(x) <= 0とh_k(x) >= 0として扱われます。[内部
        問題の解決]を参照してください。
        これはデフォルトであり、このリリースでの等式制約問題に対する唯一のオプションです。
    
        制約: 'Transform Constraints' = 'AUTO'、'NO'または
        'EQUALITIES'。
    
    'U Update Restriction' : float
        デフォルト = 0.5
    
        これは、外部反復間の(標準的な)不等式に対する
        ラグランジュ乗数の更新の境界を与える値mu_gを定義します。
        1に近い値は乗数の変化を制限し、一種の平滑化として機能し、
        低い値はより大きな変化を許容します。
    
        数値実験に基づくと、乗数の大きな変動は
        後続のステップで多数の反復を引き起こす可能性があり、
        不適切な条件設定により収束を妨げる可能性があります。
    
        特定の問題に対して、この値を実験してみる価値があるかもしれません。
        中間的な値は、より極端な値よりも推奨されます。
    
        制約: epsilon < 'U Update Restriction' < 1。
    
    'Umat Update Restriction' : float
        デフォルト = 0.3
    
        これは行列制約に対する'U Update Restriction'の等価物で、
        [概要]でmu_Aと表記されています。
        上記のアドバイスが同様に適用されます。
    
        制約: epsilon < 'Umat Update Restriction' < 1。
    
    例外
    ------
    NagValueError
        (`errno` 1)
            `handle`が初期化されていません。
    
        (`errno` 1)
            `handle`はnAG最適化モデリング
            スイートに属していないか、適切に初期化されていないか、破損しています。
    
        (`errno` 1)
            `handle`が適切に初期化されていないか、破損しています。
    
        (`errno` 2)
            このソルバーは`handle`で定義されたモデルをサポートしていません。
    
        (`errno` 3)
            問題はすでに解決中です。
    
        (`errno` 4)
            入力時、nvar = *<value>*、期待値 = *<value>*。
    
            制約: nvarは`handle`内のモデルの現在の変数数と
            一致する必要があります。
    
        (`errno` 5)
            入力時、nnzu = *<value>*。
    
            nnzuは(標準的な)制約に対するラグランジュ乗数の
            サイズと一致しません。
    
            nnzu = 0または*<value>*。
    
        (`errno` 5)
            入力時、nnzu = *<value>*。
    
            nnzuは(標準的な)制約に対するラグランジュ乗数の
            サイズと一致しません。
    
            (標準的な)制約がない場合、nnzu = 0。
    
        (`errno` 5)
            入力時、nnzua = *<value>*。
    
            nnzuaは行列制約に対するラグランジュ乗数の
            サイズと一致しません。
    
            nnzua = 0または*<value>*。
    
        (`errno` 5)
            入力時、nnzua = *<value>*。
    
            nnzuaは行列制約に対するラグランジュ乗数の
            サイズと一致しません。
    
           nnzua = 0 の場合、行列制約がありません。
    
        (`errno` 5)
            入力時、nnzuc = *<value>*。
    
            nnzucが二次錐制約のラグランジュ乗数のサイズと一致しません。
    
            二次錐制約がない場合、nnzuc = 0 となります。
    
        (`errno` 21)
            現在の開始点が使用できません。
    
        (`errno` 51)
            前処理中に問題が実行不可能であることが判明しました。
    
        (`errno` 52)
            前処理中に問題が無制限であることが判明しました。
    
        (`errno` 53)
            問題が実行不可能であると思われるため、アルゴリズムが
            停止しました。
    
        (`errno` 54)
            問題が無制限であると思われるため、アルゴリズムが
            停止しました。
    
    警告
    -----
    NagAlgorithmicWarning
        (`errno` 50)
            アルゴリズムが準最適解に収束しました。
    
            完全な精度は達成されませんでした。解は
            まだ使用可能であるはずです。
    
    NagAlgorithmicMajorWarning
        (`errno` 20)
            モニタリングステップ中にユーザーが終了を要求しました。
    
        (`errno` 22)
            外部反復の制限に達しました。
    
            要求された精度は達成されていません。
    
        (`errno` 23)
            内部サブ問題を要求された精度で解くことができませんでした。
    
            内部反復の制限に達しました。
    
        (`errno` 23)
            内部サブ問題を要求された精度で解くことができませんでした。
    
            内部サブ問題での限定的な進展が停止をトリガーしました
            (ヒューリスティックな内部停止基準)。
    
        (`errno` 23)
            内部サブ問題を要求された精度で解くことができませんでした。
    
            ライン検索または他の内部コンポーネントが失敗しました。
    
        (`errno` 24)
            進展できず、アルゴリズムが停止しました。
    
    注意
    -----
    ``handle_solve_pennon``は、ハンドルとして保存された互換性のある問題の
    ソルバーとして機能します。
    ハンドルは、問題を定義し、nAG最適化モデリングスイートの関数の
    通信手段として機能する内部データ構造を指します。
    まず、``handle_init``を呼び出してプロブレムハンドルを初期化します。
    次に、``handle_set_linobj``、
    ``handle_set_quadobj``、``handle_set_simplebounds``、
    ``handle_set_linconstr``、``handle_set_linmatineq``または
    ``handle_set_quadmatineq``のいくつかの関数を使用して、目的
    関数、(標準的な)制約、および問題の行列制約を定式化することができます。
    問題が完全に設定されたら、ハンドルをソルバーに渡すことができます。
    ハンドルが不要になったら、``handle_free``を呼び出して
    破棄し、保持されているメモリを解放する必要があります。
    nAG最適化モデリングスイートの詳細については、[E04の紹介]を参照してください。
    
    このように定義できる問題の例として、(一般的に
    非凸の)二次計画法(QP)があります
    
        [表省略]
    
    線形半正定値計画問題(SDP)
    
        [表省略]
    
   または双線形行列不等式(BMI-SDP)を含む半正定値計画問題

        [表省略]

    ここで、c、l_x、u_xはn次元ベクトル、Hは対称n*n行列、l_B、u_Bはm_B次元ベクトル、Bは一般的なm_B*n長方形行列、A_i^k、Q_ij^kは対称行列です。
    表現S⪰0は対称行列Sの固有値に対する制約を表し、すべての固有値が非負、つまり行列が半正定値であることを意味します。
    問題の定式化の詳細については、スイートの関連する関数を参照してください。

    ソルバーは、標準的な行列ペナルティ関数の適切な選択による一般化された拡張ラグランジュ法に基づいています。
    アルゴリズムの詳細な説明については、[アルゴリズムの詳細]を参照してください。
    問題に関する標準的な仮定(スレーター制約条件、実行可能集合上の目的関数の有界性、詳細はStingl(2006)を参照)の下で、アルゴリズムは局所解に収束します。
    線形SDPや凸QPなどの凸問題の場合、これは大域的解となります。
    このソルバーは、小規模な密な問題と大規模なスパースな問題の両方に適しています。

    アルゴリズムの動作とソルバーの戦略は、様々なオプション([その他のパラメータ]を参照)によって変更できます。これらは、``handle_init``によるハンドルの初期化とソルバーの呼び出しの間のいつでも、``handle_opt_set``と``handle_opt_set_file``によって設定できます。
    ソルバーが終了すると、次の解決のためにオプションを変更することができます。
    ソルバーは、様々な開始点やオプションで繰り返し呼び出すことができます。

    デフォルトの選択が'AUTO'である複数の選択肢を持つオプションがいくつかあります(例えば、'ヘッシアン密度')。
    この値は、問題の構造に基づいてソルバーにオプションの決定を委ねることを意味します。
    オプションゲッター``handle_opt_get``を呼び出して、これらのオプションの選択や他のオプションを取得することができます。

    'タスク'オプションを使用して、問題を最大化に切り替えたり、目的関数を無視して実行可能点のみを見つけたりすることができます。

    'モニター頻度'オプションを使用して、ソルバーのモニターモードをオンにすることができます。
    このモードで呼び出されたソルバーは、最適点が見つかる前でも定期的に一時停止し、呼び出しプログラムから進捗状況をモニタリングできるようにします。
    すべての重要なエラー測定値と統計が呼び出しプログラムで利用可能であり、必要に応じてソルバーを早期に終了させることができます(引数`inform`を参照)。

    ラグランジュ乗数の構造
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    アルゴリズムは内部的に、決定変数(xで表される)と標準および行列制約のラグランジュ乗数(双対変数、それぞれuとUで表される)の両方の推定値を扱います。
    初期推定値を提供したり、実行中(モニターモードがオン)に近似値を要求したり、最終値を取得したりすることができます。
    ラグランジュ乗数は2つの配列に分割され、(標準)制約の乗数uは配列`u`に、行列制約の乗数Uは配列`ua`に格納されます。
    両方の配列は制約の構造に適合する必要があります。
    
    単純な境界が定義されている場合(`handle_set_simplebounds`が正常に呼び出された場合)、`u`の最初の2n要素は対応するラグランジュ乗数に属し、各x_iの下限と上限の乗数が交互に配置されます。
    境界のいずれかが無限大に設定されている場合、対応するラグランジュ乗数は0に設定され、無視される場合があります。
    
    同様に、`u`の次の2m_B要素は、`handle_set_linconstr`によって定式化された場合、線形制約の乗数に属します。
    構成は同じで、つまり各制約の下限と上限の乗数が交互に配置され、欠落している(無限大の境界)制約にはゼロが使用されます。
    
    次元d*dの行列制約(1ブロック)のラグランジュ乗数は、同じ次元の密対称行列です。
    すべての乗数Uは、`handle_set_linmatineq`と`handle_set_quadmatineq`によって定義された行列制約と同じ順序で配列`ua`に隣接して格納されます。
    各行列の下三角部分は、パックされた列順序で格納されます([F07の紹介]を参照)。
    例えば、次元7、3、1、1の4つの行列制約がある場合、配列`ua`の次元は36になるはずです。
    最初の28要素(d_1*(d_1+1)/2)はU_1のパックされた下三角部分に属し、その後にU_2の6要素、U_3とU_4のそれぞれ1要素が続きます。
    
    ラグランジュ乗数の近似
    ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    アルゴリズムの性質上、すべての不等式ラグランジュ乗数u,Uは計算プロセス中常に正の値に保たれます。
    これは、解における非アクティブな制約のラグランジュ乗数にも適用されます。
    それらは通常厳密にゼロになるはずですが、ゼロに近い値になるだけです。
    これは、アクティブセット法に基づくソルバー(`qpconvex2_sparse_solve`など)と、`handle_solve_pennon`や内点法などの他のソルバーの結果の主な違いの1つです。
    結果として、乗数の初期推定値(提供された場合)はソルバーによって十分に正の値に調整される可能性があり、また中間の終了時に返される推定値は、すべてのカルーシュ・クーン・タッカー(KKT)条件を満たしていないため、最終値の非常に粗い近似にすぎない可能性があります。
    
    もう1つの違いは、`qpconvex2_sparse_solve`が下限と上限の不等式の乗数を1つの要素にマージし、その符号が不等式を決定することです。これは、最大で1つのアクティブな制約しかなく、非アクティブな制約の乗数は厳密にゼロになるためです。
    負の乗数は上限に関連付けられ、正の乗数は下限に関連付けられます。
    一方、`handle_solve_pennon`は両方の乗数を同時に扱うため、下限用と上限用の2つの要素で返されます([ラグランジュ乗数の構造]を参照)。
    上限の乗数を下限の乗数から引くことで、同等の結果を得ることができます。
    これは、等式が2つの不等式として解釈される場合でも成り立ちます(オプション'Transform Constraints'を参照)。
    
    参考文献
    ----------
    Ben--Tal, A and Zibulevsky, M, 1997, `Penalty/barrier multiplier
    methods for convex programming problems`, SIAM Journal on
    Optimization (7), 347--366
    
    Fujisawa, K, Kojima, M, Nakata, K, 1997, `Exploiting sparsity in
    primal-dual interior-point method for semidefinite programming`,
    Math. Programming (79), 235--253
    
    Hogg, J D and Scott, J A, 2011, `HSL MA97: a bit-compatible
    multifrontal code for sparse symmetric systems`, RAL Technical
    Report. RAL-TR-2011-024
    
    Kočvara, M and Stingl, M, 2003, `PENNON -- a code for convex
    nonlinear and semidefinite programming`, Optimization Methods and
    Software (18(3)), 317--333
    
    Kočvara, M and Stingl, M, 2007, `On the solution of large-scale SDP
    problems by the modified barrier method using iterative solvers`,
    Math. Programming (Series B) (109(2--3)), 413--444
    
    Mittelmann, H D, 2003, `An independent benchmarking of SDP and SOCP
    solvers`, Math. Programming (95), 407--430
    
    Stingl, M, 2006, `On the Solution of Nonlinear Semidefinite Programs
    by Augmented Lagrangian Methods, PhD thesis`, Institute of Applied
    Mathematics II, Friedrich--Alexander University of
    Erlangen--Nuremberg

    関連項目
    --------
    :meth:`naginterfaces.library.examples.opt.handle_solve_pennon_bmi_ex.main`
    :meth:`naginterfaces.library.examples.opt.handle_solve_pennon_lmi_ex.main`
関連情報
MENU
© 日本ニューメリカルアルゴリズムズグループ株式会社 2025
Privacy Policy  /  Trademarks