应用脊索分解提高优化性能

导言

在具有正半定约束的数学优化问题中，矩阵变量的维数往往成为限制计算效率的关键因素。然而，许多实际任务的特征在于矩阵的显着稀疏性，这开辟了使用专门分解方法的可能性。

Chordal decomposition是一种将一个大PSD约束分解为一组较小的PSD约束和线性相等形式的几个约束的方法。

*PSD约束（Positive Semidefinite constraint）是数学优化中的一种约束，要求变量的矩阵是正半定的，换句话说，矩阵描述了一种形状，如果以二次形式使用，则从不给出负值。如果原始PSD约束是稀疏的，那么分解的问题可以比原始问题更快地解决。

PSD约束是保证矩阵给出的二次形状总是非负的。这是许多系统的物理可行性的基本要求，从财务模型到工程结构。

此示例演示如何使用库。:

JuMP, MathOptChordalDecomposition, SCS, SparseArrays

以提高具有PSD限制的模型的性能。

初始数据

我们将导入并附加必要的库。

import Pkg
Pkg.add(["JuMP", "MathOptChordalDecomposition", "SCS", "SparseArrays"])
using JuMP, MathOptChordalDecomposition, SCS, SparseArrays

为了演示脊索分解的优势，我们使用来自数据的预先准备的模型。dat-s文件。

модель = read_from_file(joinpath(@__DIR__, "data.dat-s"))

A JuMP Model
├ solver: none
├ objective_sense: MIN_SENSE
│ └ objective_function_type: AffExpr
├ num_variables: 124
├ num_constraints: 1
│ └ Vector{AffExpr} in MOI.PositiveSemidefiniteConeTriangle: 1
└ Names registered in the model: none

该模型有124个决策变量和一个PSD约束。这个PSD约束是稀疏的，这意味着矩阵的许多元素都是零。

要查看矩阵，请使用函数 all_constraints 要获取限制列表，然后使用函数 constraint_object 以获得函数和集合形式的约束形式:

S = MOI.PositiveSemidefiniteConeTriangle
约束=all_constraints(model,Vector{AffExpr},S)
ogre=constraint_object（约束[1]）;
食人魔套装

MathOptInterface.PositiveSemidefiniteConeTriangle(124)

约束条件表示为向量。

食人魔[医]func

7750-element Vector{AffExpr}:
 _[1] - 0.25
 0
 _[2] - 0.5
 0
 0
 _[3] - 0.25
 0
 0
 0
 _[4]
 0
 0
 0
 ⋮
 0
 0
 0
 0
 0
 0
 0
 0
 0
 0
 0
 _[124] - 1

使用函数 reshape_vector，我们将约束转换为矩阵的形式。

F=reshape_vector（ogr。func，SymmetricMatrixShape(ogr.设置。side_dimension))

124×124 LinearAlgebra.Symmetric{AffExpr, Matrix{AffExpr}}:
 _[1] - 0.25  0           0            0     …  0           0
 0            _[2] - 0.5  0            0        0           0
 0            0           _[3] - 0.25  0        0           0
 0            0           0            _[4]     0           0
 0            0           0            0        0           0
 0            0           0            0     …  0           0
 0            0           0            0        0           0
 0            0           0.25         0        0           0
 0            0           0            0        0           0
 0            0           0            0        0           0
 0            0           0            0     …  0           0
 0            0           0            0        0           0
 0            0           0            0        0           0
 ⋮                                           ⋱              
 0            0           0            0        0           0
 0            0           0            0        0           0
 0            0           0            0        0           0
 0            0           0            0     …  0           0
 0            0           0            0        0           0
 0            0           0            0        0           0
 0            0           0            0        0           0
 0            0           0            0        0           0
 0            0           0            0     …  0           0
 0            0           0            0        0           0
 0            0           0            0        _[123] - 1  0
 0            0           0            0        0           _[124] - 1

使用函数 SparseArrays.sparse，我们变换矩阵 F 成稀疏矩阵 A:

A = SparseArrays.sparse(F)

124×124 SparseMatrixCSC{AffExpr, Int64} with 422 stored entries:
⎡⠑⢄⠂⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠂⠀⠀⡀⠀⠀⠀⠂⡀⢀⠀⠀⠀⠈⠀⠀⠀⠀⠀⠀⢀⠀⢀⠀⠀⠀⎤
⎢⠈⠀⡑⢌⠀⠠⠀⠀⠈⢀⠀⠠⠃⠀⡀⠀⠀⠀⠀⠀⠄⠀⠀⠀⠐⠠⠀⠄⠀⠀⠀⠀⠀⠠⠀⠀⠀⠀⠄⠀⎥
⎢⠀⠀⠀⡀⠑⢄⠀⠀⠀⠠⠀⠀⠀⠄⠀⠠⢈⠀⠀⠀⠠⠀⠁⠀⠄⠀⠀⠀⠀⠀⠀⠀⠀⠄⠀⠀⠀⠀⠠⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠑⢄⠀⠀⠀⠀⠀⠐⡀⠁⠀⡂⠀⠀⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠂⣀⠀⠀⠠⠠⠀⠑⎥
⎢⠀⠀⠂⢀⠀⡀⠀⠀⠑⢄⠀⠀⠡⠀⠂⠀⢀⠀⠀⠀⠀⠀⠂⠐⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠀⠀⢀⎥
⎢⠀⠀⠀⡀⠀⠀⠀⠀⠀⠀⠑⢄⠐⠔⠀⠀⠂⠁⠀⠀⠀⠂⠀⠀⠀⠀⣀⠀⠀⠀⠀⠆⠀⠂⠠⠤⠀⠆⠂⠄⎥
⎢⠀⠀⠉⠀⠀⠄⢀⠀⠁⠂⢐⠄⠑⢄⡈⠀⠄⠀⠀⠀⡄⠀⠀⠀⠈⠀⠀⠈⠀⡁⠀⠀⠀⠀⠀⠀⠙⠀⠀⠀⎥
⎢⠈⠀⠀⠈⠀⡀⠄⠈⠈⠀⠀⠀⠂⠈⠑⢄⠂⠀⠠⠀⠀⠀⠀⢄⠀⠀⠀⠀⠀⠀⢀⠀⠀⠀⠀⠈⠈⠀⠀⡀⎥
⎢⠀⠠⠀⠀⠂⠐⠠⠠⠀⠐⠌⠀⠀⠁⠈⠀⠑⢄⠀⠀⠀⠀⡀⠀⡀⠂⢀⠀⠀⠀⠀⠂⠄⠀⠀⠢⠀⠀⡀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠂⠀⠀⠑⢄⠀⠀⠀⠀⠈⠀⠀⠀⠀⠀⠀⠀⠂⠀⢀⢀⠀⠀⠈⠈⎥
⎢⠠⠀⠀⠁⠀⠂⠁⠀⠀⠀⠠⠀⠀⠉⠀⠀⠀⠀⠀⠀⠑⢄⠀⠄⠀⠠⠀⠀⠀⠁⠀⡀⠀⠈⠀⠠⠀⢀⡀⠀⎥
⎢⠀⢈⠀⠀⠁⠀⠀⠀⢈⠀⠀⠀⠀⠀⠀⢄⠀⠈⠀⠀⠀⠄⠑⢄⢀⠀⡀⠀⠀⠀⠀⠀⠈⠀⠀⠐⠀⠀⠈⠀⎥
⎢⠀⠀⠐⡀⠀⠁⠀⠀⠀⠀⠀⠀⠂⠀⠀⠀⠠⠈⠂⠀⠀⡀⠀⠐⠑⢄⠁⠐⠀⠀⠀⠀⠀⠀⠀⠈⠀⠐⠀⠀⎥
⎢⡀⠀⠀⠄⠀⠀⠀⠀⠀⠀⠀⠘⡀⠀⠀⠀⠀⠐⠀⠀⠀⠀⠀⠈⢁⠀⠑⢄⠀⠠⠐⠀⣀⠀⢠⠁⠀⢀⠂⢀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠄⠠⠀⠀⠀⠀⠀⠀⠄⠀⠀⠀⠀⠀⠀⡀⢑⢔⠀⠀⠀⠀⠀⠠⠀⠀⠀⠀⎥
⎢⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠠⠄⠀⠀⠀⠐⠠⠀⠀⠀⠀⠠⠀⠀⠀⠀⠐⠀⠀⠀⠑⢄⠀⠀⠈⠁⠀⠀⠁⡀⎥
⎢⠀⠀⠀⡀⠀⠄⠈⢠⠀⠀⠠⠀⠀⠀⠀⠀⠀⠁⠈⠀⡀⠀⠂⠀⠀⠀⠀⠘⠀⠀⠀⠀⠑⢄⠀⠀⢀⠄⠆⠀⎥
⎢⠀⠐⠀⠀⠀⠀⠀⠀⠀⠀⠀⡆⠀⠀⡀⠀⠠⡀⠀⢐⠀⡀⢀⠀⡀⠀⠄⠒⠀⡀⠆⠀⠀⠀⢑⣴⠀⠀⠀⠄⎥
⎢⠀⠐⠀⠀⠀⠀⠀⡂⠀⠐⠠⠄⠓⠀⠂⠀⠀⠀⠀⠀⠀⢀⠀⠀⢀⠀⠀⢀⠀⠀⠀⠀⠀⠔⠀⠀⠑⢄⠀⠀⎥
⎣⠀⠀⠀⠁⠀⠂⢄⠀⠀⢀⠈⠄⠀⠀⠀⠠⠀⠈⡂⠀⠀⠈⠂⠀⠀⠀⠈⢀⠀⠀⠁⠠⠈⠁⠀⠄⠀⠀⠑⢄⎦

稀疏矩阵包含422个非零元素，对应于密度为2.7%:

SparseArrays.nnz(A) / size(A, 1)^2

0.027445369406867846

计算速度

SCS是一个一阶求解器，它没有利用PSD约束的稀疏性。让我们运行这个问题的解决方案，看看花了多长时间。:

set_optimizer(model,SCS.优化器）
@时间优化！(型号)

------------------------------------------------------------------
	       SCS v3.2.10 - Splitting Conic Solver
	(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem:  variables n: 124, constraints m: 7750
cones: 	  s: psd vars: 7750, ssize: 1
settings: eps_abs: 1.0e-04, eps_rel: 1.0e-04, eps_infeas: 1.0e-07
	  alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1
	  max_iters: 100000, normalize: 1, rho_x: 1.00e-06
	  acceleration_lookback: 10, acceleration_interval: 10
	  compiled with openmp parallelization enabled
lin-sys:  sparse-direct-amd-qdldl
	  nnz(A): 124, nnz(P): 0
------------------------------------------------------------------
 iter | pri res | dua res |   gap   |   obj   |  scale  | time (s)
------------------------------------------------------------------
     0| 9.71e-01  3.87e+00  3.12e+02  2.65e+02  1.00e-01  9.67e-02 
   250| 1.20e-03  3.27e-04  2.64e-01  1.42e+02  2.92e-02  1.16e+00 
   500| 2.97e-04  7.52e-05  4.58e-02  1.42e+02  2.92e-02  2.20e+00 
   750| 2.40e-04  6.14e-05  3.89e-02  1.42e+02  2.92e-02  3.37e+00 
  1000| 1.91e-04  5.10e-05  3.30e-02  1.42e+02  2.92e-02  4.36e+00 
  1250| 1.55e-04  4.34e-05  2.78e-02  1.42e+02  2.92e-02  5.43e+00 
  1500| 1.34e-04  3.65e-05  2.34e-02  1.42e+02  2.92e-02  6.46e+00 
  1750| 2.19e-01  2.84e-01  3.74e+00  1.44e+02  2.92e-02  7.56e+00 
  2000| 1.01e-04  2.57e-05  1.67e-02  1.42e+02  2.92e-02  8.60e+00 
  2250| 9.41e-05  2.15e-05  1.42e-02  1.42e+02  2.92e-02  9.67e+00 
------------------------------------------------------------------
status:  solved
timings: total: 9.67e+00s = setup: 1.89e-02s + solve: 9.65e+00s
	 lin-sys: 4.98e-01s, cones: 8.36e+00s, accel: 3.76e-02s
------------------------------------------------------------------
objective = 141.972713
------------------------------------------------------------------
 11.288593 seconds (1.50 M allocations: 75.277 MiB, 13.45% compilation time)

我们使用SCS。优化器作为MathOptChordalDecomposition的参数。优化器。让我们启动解决方案并测量解决时间。

set_optimizer(model,()->MathOptChordalDecomposition。优化器(SCS.优化器））
@时间优化！(型号)

------------------------------------------------------------------
	       SCS v3.2.10 - Splitting Conic Solver
	(c) Brendan O'Donoghue, Stanford University, 2012
------------------------------------------------------------------
problem:  variables n: 1155, constraints m: 8781
cones: 	  z: primal zero / dual free vars: 7750
	  s: psd vars: 1031, ssize: 115
settings: eps_abs: 1.0e-04, eps_rel: 1.0e-04, eps_infeas: 1.0e-07
	  alpha: 1.50, scale: 1.00e-01, adaptive_scale: 1
	  max_iters: 100000, normalize: 1, rho_x: 1.00e-06
	  acceleration_lookback: 10, acceleration_interval: 10
	  compiled with openmp parallelization enabled
lin-sys:  sparse-direct-amd-qdldl
	  nnz(A): 2186, nnz(P): 0
------------------------------------------------------------------
 iter | pri res | dua res |   gap   |   obj   |  scale  | time (s)
------------------------------------------------------------------
     0| 2.75e+01  1.00e+00  9.93e+03 -4.81e+03  1.00e-01  1.24e-02 
   250| 6.74e-03  1.47e-03  6.73e-03  1.42e+02  1.73e+00  9.51e-01 
   350| 1.23e-04  7.46e-05  4.30e-05  1.42e+02  1.73e+00  1.34e+00 
------------------------------------------------------------------
status:  solved
timings: total: 1.34e+00s = setup: 7.11e-03s + solve: 1.33e+00s
	 lin-sys: 2.87e-01s, cones: 9.09e-01s, accel: 5.99e-03s
------------------------------------------------------------------
objective = 141.990174
------------------------------------------------------------------
  7.378757 seconds (8.81 M allocations: 443.399 MiB, 1.60% gc time, 81.60% compilation time)

现在决定时间不到一秒钟。这种性能差异是由于使用了脊索分解。分解的任务引入了新的变量（现在有1155个而不是124个）和约束（现在有115个PSD约束而不是一个），但每个PSD约束都变得比原来的小。

分解=不安全后端（模型）

MathOptChordalDecomposition.Optimizer{CliqueTrees.MF}
├ ObjectiveSense: MIN_SENSE
├ ObjectiveFunctionType: MOI.ScalarAffineFunction{Float64}
├ NumberOfVariables: 1155
└ NumberOfConstraints: 116
  ├ MOI.VectorAffineFunction{Float64} in MOI.Zeros: 1
  └ MOI.VectorOfVariables in MOI.PositiveSemidefiniteConeTriangle: 115

让我们计算每个尺寸的PSD约束的数量:

数量=Dict{Int,Int}()
为了moi的ci。得到（分解，MOI。ListOfConstraintIndices{MOI.VectorOfVariables，S}（））
    set=MOI。得到（分解，MOI。约束集(),ci)
    n = set.side_dimension
    数量[n]=得到(数量,n,0)+1
end
数量

Dict{Int64, Int64} with 10 entries:
  5  => 7
  4  => 15
  6  => 3
  7  => 3
  2  => 33
  10 => 2
  9  => 2
  8  => 3
  3  => 35
  1  => 12

最大的PSD约束现在的大小为10，显着小于原始的124x124矩阵。

结论

该研究证明了脊索分解作为一种加速稀疏矩阵约束半定规划问题求解的方法的有效性。

进一步研究的前景包括使该方法适应混合结构（PSD约束与其他类型约束的组合）的问题，优化为特殊类稀疏矩阵构建脊索扩展的算法，以及根据对问题结构的分析制定评估分解可行性的标准。

该方法的实施确保了结果的可重复性，并为将脊索分解集成到处理大规模优化任务的研究人员和工程师的工作流程中创造了基础。