Processing math: 100%

半正定値計画法(SDP)を用いた行列補完

nAG Library for Python Example集

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

半正定値計画法(SDP)を用いた行列補完

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

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

はじめに

nr人の回答者とnq個の質問からなる調査から始めます。残念ながら、いくつかのエントリが欠落しています…

以下はnr=15nq=6の例です ˆY=[0.40.60.40.81.00.80.20.80.21.00.40.00.20.40.20.21.00.80.20.61.00.21.00.40.60.01.00.41.00.20.20.40.41.01.00.81.00.20.61.00.20.60.20.4] 欠損値を埋める最良の方法は何でしょうか?一つの見方として、類似した回答をした回答者は類似したままであるべきだと言えるでしょう。実際には、これはYのランクを最小化することでモデル化できます。 minYrank(Y)Yij=ˆYij ランク最小化は一般的にNP困難ですが、ヒューリスティックによって近似することができます。それは行列の核ノルムを最小化することです。行列の核ノルムは、その特異値の和です。ランク不足の行列は(複数の)ゼロ特異値を持つ必要があります。特異値は常に非負であるという事実を考えると、核ノルムの最小化は圧縮センシングにおけるl1ノルムと同じ効果を持ちます。つまり、多くの特異値をゼロにすることを促進し、したがって元のランク最小化問題のヒューリスティックとみなすことができます。 minY||Y||Yij=ˆYij これは以下のSDP問題を解くことで達成できます: minY,W1,W2trace(W1)+trace(W2)Yij=ˆYij[W1YYTW2]0

# データを定義することから始めます
nr = 15
nq = 6
# 行列、欠損値は-1に設定
Y = np.array([[-1.0,-1.0,-1.0,-1.0,-1.0, 0.4],
                [ 0.6, 0.4, 0.8,-1.0,-1.0,-1.0],
                [-1.0,-1.0, 0.8,-1.0, 0.2,-1.0],
                [ 0.8, 0.2,-1.0,-1.0,-1.0,-1.0],
                [-1.0, 0.4,-1.0, 0.0,-1.0, 0.2],
                [ 0.4,-1.0,-1.0, 0.2,-1.0, 0.2],
                [-1.0, 0.8, 0.2, 0.6,-1.0,-1.0],
                [-1.0,-1.0, 0.2,-1.0,-1.0,-1.0],
                [-1.0, 0.4,-1.0, 0.6, 0.0,-1.0],
                [-1.0,-1.0, 0.4,-1.0,-1.0,-1.0],
                [-1.0,-1.0, 0.2, 0.2, 0.4, 0.4],
                [-1.0,-1.0,-1.0,-1.0, 1.0, 0.8],
                [ 1.0,-1.0, 0.2,-1.0,-1.0, 0.6],
                [-1.0,-1.0,-1.0,-1.0,-1.0, 0.2],
                [ 0.6,-1.0, 0.2, 0.4,-1.0,-1.0]])

# 欠落要素の数を数える
n_miss = 0
for i in range(nr):
    for j in range(nq):
        if Y[i][j] <= -1.0:
            n_miss += 1

これでSDPの問題定義を始めることができます。W1W2はそれぞれnr×nrnq×nqの完全な変数行列です。さらに、欠損値も最適化問題の変数となります。

# 変数の数: 欠損値 + 2つの完全行列 nrxnr と nqxnq 
# (上三角部分のみが必要です)
nvar = int(n_miss + nr*(nr+1)/2 + nq*(nq+1)/2)

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

# 線形目的関数を0.0に初期化する
cvec = np.zeros(nvar)

線形行列不等式の定義を始めます。まず、ˆYの既知の値を、行列不等式の定数部分に格納します。

# スパース行列不等式行列を初期化する 
nnzasum = nvar + nr*nq - n_miss
nnza = np.empty(nvar+1, dtype=int)
irowa = np.empty(nnzasum, dtype=int)
icola = np.empty(nnzasum, dtype=int)
a = np.empty(nnzasum, dtype=float)

# 定数ブロック内のゼロでない要素の数
nnza[0] = nr*nq - n_miss

# 組み合わせの中の他のすべての行列は1つの要素しか持ちません
nnza[1:] = 1

# Y の既知の値で A_0 を埋める
idx = 0
for i in range(nr):
    for j in range(nq):
        if Y[i][j] >= 0.0:
            irowa[idx] = i + 1
            icola[idx] = nr + j + 1
            a[idx] = -Y[i][j]
            idx += 1
次に上左のブロックを埋めます [W1]

これはW1の上三角部分で構成されています。

# 1x1のマトリックスブロック
# 目的ベクトルを同時に埋める
idxobj = 0
for i in range(nr):
    # これはW1の対角要素で、トレースに属します(目的を設定)
    cvec[idxobj] = 1.0
    for j in range(i, nr):
        irowa[idx] = i + 1
        icola[idx] = j + 1
        a[idx] = 1.0
        idx += 1
        idxobj += 1
右下のブロックの時間 [W2]

W2の上三角部分から成る。

for i in range(nq):
    # これはW2の対角要素で、トレースに属します(目的を設定)
    cvec[idxobj] = 1.0
    for j in range(i, nq):
        irowa[idx] = nr + i + 1
        icola[idx] = nr + j + 1
        a[idx] = 1.0
        idx += 1
        idxobj += 1

最後に、Yの欠損要素

for i in range(nr):
    for j in range(nq):
        if Y[i][j] < 0.0:
            irowa[idx] = i + 1 
            icola[idx] = nr + j + 1
            a[idx] = 1.0
            idx +=1

問題のハンドルを更新する

# ハンドルに線形目的関数を追加する
opt.handle_set_linobj(handle, cvec)

# 線形行列不等式を追加する
dima = nr + nq
idblk = opt.handle_set_linmatineq(handle, dima, nnza, irowa, icola, a, blksizea=None, idblk=0)

# オプション引数を設定する
for option in ['Print Options = No',
               'Initial X = Automatic',
               'Dimacs = Check']:
    opt.handle_opt_set(handle, option)

問題を解く準備ができました

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

# ソルバーを呼び出す
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                196
     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              21
     detected normal inq.           0
       linear                       0
       nonlinear                    0
   Objective function          Linear

 --------------------------------------------------------------
  it|  objective |  optim  |   feas  |  compl  | pen min |inner
 --------------------------------------------------------------
   0  0.00000E+00  9.17E+01  1.91E+00  0.00E+00  2.00E+00   0
   1  1.53207E+02  6.26E-03  0.00E+00  1.50E+02  2.00E+00  13
   2  7.10956E+01  1.75E-02  0.00E+00  6.45E+01  9.04E-01   4
   3  3.68500E+01  1.36E-02  0.00E+00  2.73E+01  4.08E-01   5
   4  2.25150E+01  1.07E-02  0.00E+00  1.15E+01  1.85E-01   4
   5  1.64698E+01  1.10E-02  0.00E+00  4.84E+00  8.34E-02   5
   6  1.38820E+01  4.07E-03  0.00E+00  2.06E+00  3.77E-02   5
   7  1.27611E+01  7.27E-03  0.00E+00  8.83E-01  1.70E-02   5
   8  1.22745E+01  4.16E-03  0.00E+00  3.82E-01  7.70E-03   4
   9  1.20628E+01  3.53E-03  0.00E+00  1.65E-01  3.48E-03   4
  10  1.19705E+01  3.14E-03  0.00E+00  7.20E-02  1.57E-03   4
  11  1.19303E+01  3.67E-03  0.00E+00  3.13E-02  7.11E-04   4
  12  1.19127E+01  3.65E-03  0.00E+00  1.37E-02  3.21E-04   4
  13  1.19050E+01  4.02E-03  0.00E+00  5.96E-03  1.45E-04   4
  14  1.19017E+01  1.35E-04  0.00E+00  2.60E-03  6.56E-05   5
 --------------------------------------------------------------
  it|  objective |  optim  |   feas  |  compl  | pen min |inner
 --------------------------------------------------------------
  15  1.19002E+01  9.38E-07  0.00E+00  1.13E-03  2.96E-05   5
  16  1.18996E+01  1.51E-06  0.00E+00  4.90E-04  1.34E-05   5
  17  1.18993E+01  1.45E-10  0.00E+00  2.12E-04  6.06E-06   6
  18  1.18992E+01  3.23E-10  0.00E+00  9.14E-05  2.74E-06   6
  19  1.18991E+01  1.09E-09  0.00E+00  3.93E-05  1.24E-06   6
  20  1.18991E+01  1.70E-09  0.00E+00  1.68E-05  5.59E-07   6
  21  1.18991E+01  6.16E-09  0.00E+00  7.15E-06  2.53E-07   6
  22  1.18991E+01  1.55E-08  0.00E+00  3.02E-06  1.14E-07   6
  23  1.18991E+01  2.19E-08  0.00E+00  1.27E-06  5.16E-08   6
 --------------------------------------------------------------
 Status: converged, an optimal solution found
 --------------------------------------------------------------
 Final objective value                1.189910E+01
 Relative precision                   1.359557E-07
 Optimality                           2.187594E-08
 Feasibility                          0.000000E+00
 Complementarity                      1.267053E-06
 DIMACS error 1                       1.093797E-08
 DIMACS error 2                       0.000000E+00
 DIMACS error 3                       0.000000E+00
 DIMACS error 4                       0.000000E+00
 DIMACS error 5                       5.004717E-08
 DIMACS error 6                       5.109456E-08
 Iteration counts
   Outer iterations                             23
   Inner iterations                            122
   Linesearch steps                            332
 Evaluation counts
   Augm. Lagr. values                          146
   Augm. Lagr. gradient                        146
   Augm. Lagr. hessian                         122
 --------------------------------------------------------------

これで完成した行列を取得できます!

# 行列の欠けている要素を埋めてください
idx = int(nr*(nr+1)/2 + nq*(nq+1)/2)
for i in range(nr):
    for j in range(nq):
        if Y[i][j] < 0.0:
            Y[i][j] = x[idx]
            idx += 1
            
np.set_printoptions(precision=1)
print(Y)

rk = np.linalg.matrix_rank(Y, tol=1.0e-05)
print('The rank of the filled matrix is:', rk)
[[0.5 0.3 0.2 0.2 0.4 0.4]
 [0.6 0.4 0.8 0.2 0.3 0.4]
 [0.4 0.3 0.8 0.  0.2 0.2]
 [0.8 0.2 0.3 0.4 0.3 0.4]
 [0.  0.4 0.2 0.  0.2 0.2]
 [0.4 0.1 0.2 0.2 0.1 0.2]
 [0.6 0.8 0.2 0.6 0.2 0.4]
 [0.1 0.1 0.2 0.  0.  0.1]
 [0.6 0.4 0.1 0.6 0.  0.3]
 [0.2 0.1 0.4 0.  0.1 0.1]
 [0.5 0.3 0.2 0.2 0.4 0.4]
 [0.7 0.4 0.3 0.  1.  0.8]
 [1.  0.3 0.2 0.5 0.5 0.6]
 [0.2 0.1 0.1 0.1 0.2 0.2]
 [0.6 0.3 0.2 0.4 0.2 0.3]]
埋められた行列のランクは: 4    
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'を参照)。

    References
    ----------
    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