単精度の quadprog コードの生成

この例では、quadprog ソルバーを使用して単精度ターゲットのコードを生成する方法を示します。二次計画問題は、$$L^1$$ ペナルティを使用して二次問題の比較的スパースな解を生成することです。この手法は LASSO 回帰と呼ばれ、lasso (Statistics and Machine Learning Toolbox)関数で使用できます。

具体的には、与えられた行列 $$A$$、ベクトル $$b$$、およびペナルティ パラメーター $$\lambda$$ について、LASSO 回帰は次の問題を解きます。

$$\underset{x}{\min }\frac{1}{2}\Vert Ax-b{\Vert }_2^2+\lambda \Vert x{\Vert }_1$$.

$$L^1$$ ペナルティのため、LASSO 回帰の結果には通常、ゼロのエントリがいくつか含まれます。

この問題を、quadprog による解法に適した二次計画法の形式に再キャストします。$$e$$ は、長さ n の 1 の列ベクトルを表すとします。n は $$A$$ の列数です。$$s$$ は長さ n の列ベクトルを表すとすると、$$e^{T}s$$ は $$e$$ と $$s$$ のドット積です。次の二次計画法は、元の LASSO 問題と同等です。

$$\begin{array}{l}\underset{x,z,s}{\min }\frac{1}{2}{z}^{T}z+\lambda e^{T}s\\ \text{subject}\text{to}\\ Ax-b=z\\ -s\le x\le s\\ s\ge 0.\end{array}$$

$$A$$ 係数、$$b$$ 係数、および $$\lambda$$ 係数の単精度データを作成します。$$A$$ を 30 行 8 列の疑似乱数行列として設定します。$$b$$ を、エントリの約半分がゼロに等しいベクトル x_exp から作成された 30 行 1 列の疑似乱数ベクトルとして設定します。$$\lambda$$ を単精度スカラー 1 として作成します。

rng default;
m = 30;
n = 8;

datatype = "single";
x_exp = randi([0 1],n,1,datatype) .* randn(n,1,datatype);
A = rand(m,n,datatype);
b = A*x_exp;
lambda = cast(1,datatype);

問題を単精度で解く関数の作成

以下に示す solveLassoProb 関数コードをファイルとして MATLAB® パス上に保存します。

function [x,reserror,flag] = solveLassoProb(A,b,lambda)

[m,n] = size(A);

% Build the objective.
H = blkdiag(zeros(n,"like",b), ...  % x
            eye(m,"like",b), ...    % z
            zeros(n,"like",b));     % s

f = [zeros(n,1,"like",b);           % x
     zeros(m,1,"like",b);           % z
     lambda*ones(n,1,"like",b)];    % s

% Define z as the residual (Ax - b).
Aeq = [A  -eye(m,m,"like",b) zeros(m,n,"like",b)];
beq = b;

% Add inequalities between x and s. Ensure datatype consistency.
Aineq = [eye(n,"like",b) zeros(n,m,"like",b) -eye(n,"like",b)
         -eye(n,"like",b) zeros(n,m,"like",b) -eye(n,"like",b)];
bineq = zeros(2*n,1,"like",b);

% s must be nonnegative.
lb = [-optim.coder.infbound(n+m,1,"like",b); zeros(n,1,"like",b)];
ub = optim.coder.infbound(2*n+m,1,"like",b);

x0 = zeros(2*n+m,1,"like",b);

% Set options to use the "active-set" algorithm (required for code
% generation) and UseCodegenSolver=true to have the resulting code be the
% same whether run in MATLAB or by using code generation.
options = optimoptions("quadprog",Algorithm="active-set",UseCodegenSolver=true);

[xall,~,flag] = quadprog(H,f,Aineq,bineq,Aeq,beq,lb,ub,x0,options);
x = xall(1:n);
reserror = norm(A*x - b,2);
end

問題を倍精度で解く

単精度コード生成を使用して問題を解く前に、倍精度係数を使用して問題を解くことで、問題の定式化が正しいことを確認します。

Ad = double(A);
bd = double(b);
lambdad = double(lambda);
[x,reserror,flag] = solveLassoProb(Ad,bd,lambdad);

Minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

解の変数 x と元の変数 x_exp を比較します。

array2table([x,x_exp],"VariableNames",["LASSO Solution","Original Variables"])

ans=8×2 table
    LASSO Solution    Original Variables
    ______________    __________________

          3.5587             3.5784     
          2.7068             2.7694     
      5.0182e-16                  0     
          2.9592             3.0349     
          0.5613             0.7254     
      -1.102e-17                  0     
      1.7786e-16                  0     
     -2.3311e-16           -0.20497     

LASSO 回帰は、変数 5 と変数 8 に大きく影響しています。全体的な残差誤差はかなり大きくなります。

disp(reserror)

    0.4694

より正確な解を求めるには、ペナルティ項 $$\lambda$$ を 0.025 まで下げます。

lambdad = double(0.025);
[x,reserror,flag] = solveLassoProb(Ad,bd,lambdad);

Minimum found that satisfies the constraints.

Optimization completed because the objective function is non-decreasing in 
feasible directions, to within the value of the optimality tolerance,
and constraints are satisfied to within the value of the constraint tolerance.

array2table([x,x_exp],"VariableNames",["LASSO Solution","Original Variables"])

ans=8×2 table
    LASSO Solution    Original Variables
    ______________    __________________

          3.5765             3.5784     
          2.7621             2.7694     
      7.6581e-16                  0     
          3.0298             3.0349     
         0.71407             0.7254     
      1.7349e-16                  0     
      3.7139e-16                  0     
        -0.17967           -0.20497     

変数 5 と変数 8 がターゲット (元の) 値に近くなっています。残差の値を確認します。

disp(reserror)

    0.0357

残差は以前の約 1/10 に減少しています。$$\lambda$$ の新しい値を使用して単精度コードを生成します。

単精度コードの生成と問題の求解

codegen を呼び出して MEX ファイル ターゲットを作成します。

argList = {A,b,lambda};
codegen -config:mex -args argList solveLassoProb

Code generation successful.

MATLAB により、現在のフォルダーにファイル solveLassoProb_mex が作成されます。

生成された単精度コードを使用して、問題を解きます。

lambda = single(lambdad);
[xs,reserrors,flag] = solveLassoProb_mex(A,b,lambda);

単精度の結果を倍精度の結果と元のデータの両方と比較します。

array2table([xs,x,x_exp],"VariableNames",["Codegen Single Solution","Double Solution","Original Variables"])

ans=8×3 table
    Codegen Single Solution    Double Solution    Original Variables
    _______________________    _______________    __________________

               3.5765                3.5765              3.5784     
               2.7621                2.7621              2.7694     
          -4.4278e-08            7.6581e-16                   0     
               3.0298                3.0298              3.0349     
              0.71407               0.71407              0.7254     
           1.3336e-07            1.7349e-16                   0     
           7.4389e-08            3.7139e-16                   0     
             -0.17967              -0.17967            -0.20497     

単精度の解の変数は、ゼロに近い値を除いて、表示精度に対して倍精度の解の変数と同じです。

参考

quadprog | optimwarmstart

トピック

• quadprog 用のコード生成の背景
• quadprog のコード生成