主要内容

本页采用了机器翻译。点击此处可查看最新英文版本。

具有内点算法的大型稀疏二次规划

此示例显示了遇到稀疏问题时使用稀疏算法的价值。矩阵有 n 行,其中您选择 n 作为较大值,以及一些非零对角线带。大小为 n×n 的满矩阵会用尽所有可用内存,但稀疏矩阵则不会出现问题。

问题是最小化 x'*H*x/2 + f'*x,但要满足

x(1) + x(2) + ... + x(n) <= 0,

其中 f = [-1;-2;-3;...;-n]H 是稀疏对称带状矩阵。

创建稀疏二次矩阵

根据向量 [3,6,2,14,2,6,3] 的移位创建一个对称循环矩阵,其中 14 位于主对角线上。让矩阵为 n×n,其中 n = 30,000

n = 3e4;
H2 = speye(n);
H = 3*circshift(H2,-3,2) + 6*circshift(H2,-2,2) + 2*circshift(H2,-1,2)...
    + 14*H2 + 2*circshift(H2,1,2) + 6*circshift(H2,2,2) + 3*circshift(H2,3,2);

查看矩阵结构。

spy(H)

创建线性约束和目标

线性约束是解元素之和为非正值。目标函数包含一个以向量 f 表示的线性项。

A = ones(1,n);
b = 0;
f = 1:n;
f = -f;

求解问题

使用 interior-point-convex' 算法求解二次规划问题。为了防止求解器过早停止,请将 StepTolerance 选项设置为 0

options = optimoptions(@quadprog,'Algorithm','interior-point-convex','StepTolerance',0);
[x,fval,exitflag,output,lambda] = ...
    quadprog(H,f,A,b,[],[],[],[],[],options);
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.

<stopping criteria details>

在许多计算机上,当 n = 30,000 时,您无法创建完整的 n×n 矩阵。因此,您仅使用稀疏矩阵即可运行此问题。

检查解

查看与线性不等式相关的目标函数值、迭代次数和拉格朗日乘数。

fprintf('The objective function value is %d.\nThe number of iterations is %d.\nThe Lagrange multiplier is %d.\n',...
    fval,output.iterations,lambda.ineqlin)
The objective function value is -3.133073e+10.
The number of iterations is 7.
The Lagrange multiplier is 1.500050e+04.

由于没有下界、上界或线性等式约束,唯一有意义的拉格朗日乘数是 lambda.ineqlin。因为 lambda.ineqlin 非零,所以您可以知道不等式约束是有效的。评估约束以查看解是否在边界上。

fprintf('The linear inequality constraint A*x has value %d.\n',A*x)
The linear inequality constraint A*x has value 9.150244e-08.

解分量的总和为零,在公差范围内。

x 有三个区域:初始部分、最终部分和覆盖大部分解的近似线性部分。绘制三个区域。

subplot(3,1,1)
plot(x(1:60))
title('x(1) Through x(60)')
subplot(3,1,2)
plot(x(61:n-60))
title('x(61) Through x(n-60)')
subplot(3,1,3)
plot(x(n-59:n))
title('x(n-59) Through x(n)')

另请参阅

|

主题