0

私は Python を使用してダンツィヒ セレクターをプログラミングする任務を負っていますが、ガイドラインは与えられておらず、線形計画法やデータ サイエンスの経験もあまりありません。LP モジュールのマニュアルや、このサイトの他の質問で必要な情報が見つかりません。

これが問題です。列ベクトル ^β を探しています。この問題に対処する方法がわからないため、プログラムのこの部分のコードがありません。いくつかのアプローチを試しましたが、問題を正しく反映したものはなかったため、拒否して削除しました。

分||^β||l1 st ||xT(y-x^β(||l(inf) <= δ)

線形プログラムとして書き直すことができます。

^β は kx1 列ベクトルであり、私が探しているダンツィヒ セレクターです。

  • y は観測値/応答の nx1 列ベクトルです
  • X は nxk サンプル行列で、ここで k >> n
  • δ はノイズ変数です

役立つ可能性のある詳細を次に示します。

これまでの作業コードは次のとおりです。データ値はすべて単なるサンプル/プレースホルダーです。X、y、およびいくつかの δ の値を既に用意しています。しかし、^β を与える適切な LP 関数を見つけることができません。

import numpy as np
import random
import math

#n = no. runs = 5
n = 5
#k = no. variables = 23
k = 23

#y = vector of observations/responses (nx1, binary decisions)
y = np.array([[1],
              [0],
              [0],
              [1],
              [0]])
#X = predictor/sample matrix (nxk)
X = np.array([[1.1, 0, 0.7, 0.8, 0.9, 0.2, 0.3, 0.5, 0.2, 0.2, 1.2, 1.1, 0.5, 0.5, 0.7, 1.2, 1.3, 0.8, 0.9, 1.7, 1.2, 1.9, 0.9],
              [0.3, 0.1, 0.7, 0.4, 0.9, 0.9, 0.1, 0.8, 0.1, 0.2, 1.1, 0, 0.9, 0.4, 1.4, 1.4, 0.1, 0.5, 1.8, 1.6, 1.2, 1.8, 0.3],
              [0.1, 0.1, 0.3, 0.9, 0.7, 0.8, 0, 0.7, 0.8, 0.2, 1.1, 1.1, 0.5, 0.5, 0.8, 1.5, 0.2, 0.5, 1.6, 1.5, 1.2, 1.7, 0.5],
              [1.2, 0.2, 0.9, 0.8, 0.6, 0.2, 0.3, 0.5, 0.3, 0.2, 1.2, 1.1, 0.5, 0, 0.7, 1.2, 1.3, 0.8, 0.9, 1.7, 1.2, 1.9, 0.9],
              [0.2, 0.1, 0.6, 0, 0.5, 1.1, 0.2, 0.5, 0.9, 0.2, 1.2, 1.1, 0.8, 1.6, 0.5, 1.3, 0.2, 0.5, 1.7, 1.2, 1.2, 1.9, 0.1]])


#estimate missing data (0)
X_row_minima = np.where(X>0,X,X.max()).min(1)
X[X==0] = X_row_minima/2

#unit length normalize X
X = X/np.linalg.norm(X, ord=2, axis=1, keepdims=True)

#standardize y to zero mean
y = y - np.mean(y) / np.std(y)

#transpose X (kxn)
Xt = np.transpose(X)

#solve d0
Xty = np.matmul(Xt,y)
d0 = max(abs(Xty))

#generate 100 evenly-spaced d values
d = np.linspace(0, d0, 100)

これは、このサイトでの私の最初の投稿です。他の投稿と比較して、投稿の詳細が不足していることをお詫び申し上げます。

4

1 に答える 1

0

これまで Danzig Selector を理解する必要がなかったため、これが正しいかどうかはわかりません。可能であれば、回答に何が期待できるかがわかっているデータセットでテストすることをお勧めします。私が言えることから、あなたの問題はこのようなものです。

minimize:   sum(u)
subject to: -u_i <= beta_i for all i in len(u)
            u_i >= beta_i for all i in len(u)
            [X^T (y - X beta)]_j <= delta_j for all j in len(y)
            [X^T (y - X beta)]_j >= -delta_j for all j in len(y)

これを書くのcvxpyはかなり簡単だと思います。私が離れている場合、うまくいけば、それはあなたを正しい方向に導きます. の値を使用するX, y, k:

import cvxpy as cp
beta = cp.Variable((k,1))
u = cp.Variable((k,1))
delta = 1

objective = cp.Minimize(cp.sum(u))

constraints = [beta >= -u,
               beta <= u,
               Xt @ (y - X @ beta) <= delta,
               Xt @ (y - X @ beta) >= -delta]
               

prob = cp.Problem(objective, constraints)
prob.solve(verbose=True)

print("Problem status: ", prob.status)
print("Optimal value", prob.value)
print("Beta values", beta.value)

私の結果の出力:

===============================================================================
                                     CVXPY
                                    v1.1.15
===============================================================================
(CVXPY) Oct 12 11:00:50 PM: Your problem has 46 variables, 4 constraints, and 0 parameters.
(CVXPY) Oct 12 11:00:50 PM: It is compliant with the following grammars: DCP, DQCP
(CVXPY) Oct 12 11:00:50 PM: (If you need to solve this problem multiple times, but with different data, consider using parameters.)
(CVXPY) Oct 12 11:00:50 PM: CVXPY will first compile your problem; then, it will invoke a numerical solver to obtain a solution.
-------------------------------------------------------------------------------
                                  Compilation
-------------------------------------------------------------------------------
(CVXPY) Oct 12 11:00:50 PM: Compiling problem (target solver=ECOS).
(CVXPY) Oct 12 11:00:50 PM: Reduction chain: Dcp2Cone -> CvxAttr2Constr -> ConeMatrixStuffing -> ECOS
(CVXPY) Oct 12 11:00:50 PM: Applying reduction Dcp2Cone
(CVXPY) Oct 12 11:00:50 PM: Applying reduction CvxAttr2Constr
(CVXPY) Oct 12 11:00:50 PM: Applying reduction ConeMatrixStuffing
(CVXPY) Oct 12 11:00:50 PM: Applying reduction ECOS
(CVXPY) Oct 12 11:00:50 PM: Finished problem compilation (took 3.027e-03 seconds).
-------------------------------------------------------------------------------
                                Numerical solver
-------------------------------------------------------------------------------
(CVXPY) Oct 12 11:00:50 PM: Invoking solver ECOS  to obtain a solution.

ECOS 2.0.7 - (C) embotech GmbH, Zurich Switzerland, 2012-15. Web: www.embotech.com/ECOS

It     pcost       dcost      gap   pres   dres    k/t    mu     step   sigma     IR    |   BT
 0  +0.000e+00  -1.596e+02  +3e+02  2e-01  6e-01  1e+00  4e+00    ---    ---    1  1  - |  -  -
 1  +4.462e+00  -6.924e+00  +4e+01  1e-02  6e-02  2e-01  4e-01  0.9228  4e-02   0  0  0 |  0  0
 2  +4.980e-01  -8.219e-01  +4e+00  1e-03  8e-03  3e-02  5e-02  0.9010  3e-02   0  0  0 |  0  0
 3  +3.673e-02  -4.271e-02  +3e-01  9e-05  5e-04  1e-03  3e-03  0.9630  2e-02   0  0  0 |  0  0
 4  +4.054e-04  -4.749e-04  +3e-03  1e-06  5e-06  1e-05  3e-05  0.9890  1e-04   0  0  0 |  0  0
 5  +4.490e-06  -5.272e-06  +3e-05  1e-08  6e-08  2e-07  4e-07  0.9890  1e-04   0  0  0 |  0  0
 6  +4.972e-08  -5.852e-08  +4e-07  1e-10  6e-10  2e-09  4e-09  0.9890  1e-04   0  0  0 |  0  0
 7  +5.507e-10  -6.496e-10  +4e-09  1e-12  7e-12  2e-11  4e-11  0.9890  1e-04   0  0  0 |  0  0

OPTIMAL (within feastol=6.9e-12, reltol=-nan(ind), abstol=4.1e-09).
Runtime: 0.002074 seconds.

-------------------------------------------------------------------------------
                                    Summary
-------------------------------------------------------------------------------
(CVXPY) Oct 12 11:00:50 PM: Problem status: optimal
(CVXPY) Oct 12 11:00:50 PM: Optimal value: 5.507e-10
(CVXPY) Oct 12 11:00:50 PM: Compilation took 3.027e-03 seconds
(CVXPY) Oct 12 11:00:50 PM: Solver (including time spent in interface) took 2.358e-03 seconds
Problem status:  optimal
Optimal value 5.507030257315668e-10
Beta values [[ 1.43866286e-11]
 [ 2.52691994e-12]
 [ 1.26496527e-11]
 [ 1.25889902e-11]
 [ 1.32324154e-11]
 [ 9.34898130e-12]
 [ 4.19460380e-12]
 [ 1.07283774e-11]
 [ 7.54676427e-12]
 [ 3.94492714e-12]
 [ 1.68905798e-11]
 [ 1.61644206e-11]
 [ 1.09177395e-11]
 [ 9.11945636e-12]
 [ 1.38073390e-11]
 [ 1.51133406e-11]
 [ 1.59259478e-11]
 [ 1.24992988e-11]
 [ 1.10965990e-11]
 [ 1.50106084e-11]
 [ 1.66261028e-11]
 [-1.92542442e-13]
 [ 1.24499858e-11]]
于 2021-10-13T05:07:34.460 に答える