生成单精度 `quadprog` 代码

此示例使用:

此示例展示如何使用 quadprog 求解器为单精度目标生成代码。二次规划问题是利用 $L^{1}$ 惩罚来生成二次问题的相对稀疏的解。这种技术称为 Lasso 回归，可用于lasso (Statistics and Machine Learning Toolbox)功能。

具体来说，给定矩阵 $A$ 、向量 $b$ 和惩罚参数 $λ$ ，Lasso 回归可以解决问题

$\min_{x} \frac{1}{2} ‖ A x - b ‖_{2}^{2} + λ ‖ x ‖_{1}$ .

由于 $L^{1}$ 惩罚，Lasso 回归的结果通常有几个零条目。

将此问题重新转换为适合通过 quadprog 解决的二次规划形式。令 $e$ 表示长度为 n 的列向量，其中 n 是 $A$ 中的列数。令 $s$ 表示长度为 n 的列向量，因此 $e^{T} s$ 是 $e$ 和 $s$ 的点积。以下二次规划等同于原始 Lasso 问题：

$\begin{array}{l} \min_{x, z, s} \frac{1}{2} z^{T} z + λ e^{T} s \\ subject to \\ A x - b = z \\ - s \leq x \leq s \\ s \geq 0 . \end{array}$

为 $A$ 、 $b$ 和 $λ$ 系数创建单精度数据。将 $A$ 设置为 30×8 的伪随机矩阵。将 $b$ 设置为由向量 x_exp 构成的 30×1 伪随机向量，该向量大约有一半的元素等于零。将 $λ$ 创建为单精度标量 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

套索回归在很大程度上影响变量 5 和 8。残差误差总体来说比较大。

disp(reserror)

    0.4694

为了获得更准确的解，将惩罚项 $λ$ 降低至 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

残差大约比以前小了 10 倍。使用新的 $λ$ 值生成单精度代码。

生成单精度代码并解决问题

调用 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` 代码