Large-Scale SDP

1

Hello.
I am working on an SDP relaxation problem. When the size of the positive semidefinite variable is small (2828), the computation time is okay. However, when I tested a large-scale system (236236 for the positive semidefinite variable), the computation time is really high.
Is there any good solver for large-scale PSD? ( I am using mosek. SDPT3 is also not good)
Or is there any other way to model the PSD constraint (W>=0)?

I should mention that based on profile view, a lot of computation time is spent on the .mtime function.
Thanks.

2

Perhaps you can show the Mosek output. Then one of the Mosek people who read the post might be able to say something.

3

Let n be the dimension of largest matrix variable.

The computational complexity of iteration is at least n^3. Now that means n=2828 and n=236236 is too vastly different sized problems.

Mosek will spend along time on a problem n=236236 (and I am surprised you would expect otherwise).

4

@Michal_Adamaszek thinks you mean n=28 and n=236. If that is the case you are in the small scale world. So we need logs to say some something relevant.

5

Thank you for the reply.
Code:
cvx_begin sdp
cvx_solver mosek

  variable W(2*length(Nbus),2*length(Nbus)) 
  variable obj_gen(1,length(Gloc))        
  minimize sum(obj_gen)
  subject  to

%% constraint related to upper and lower bound for generated active ad reactive power

      for obj_gen_idx=1:length(Gloc)
          obj_gen1=Gloc(obj_gen_idx);
          obj_gen(obj_gen_idx)==100*trace(Y_k(obj_gen1)*W);
      end


      for g_idx = 1:length(Nbus)
          100*trace(Y_k(g_idx)*W)<=(Pmax(g_idx)-PD(g_idx)) ;
          100*trace(Y_k(g_idx)*W)>=(Pmin(g_idx)-PD(g_idx));
          100*trace(Y_k_bar(g_idx)*W)<=Qmax(g_idx)-QD(g_idx);
          100*trace(Y_k_bar(g_idx)*W)>=Qmin(g_idx)-QD(g_idx);
      end

%% semidefinite Constraint

W>=0;

cvx_end

where length(Nbus)=118 and length(Gloc)=54.

Mosek log:
Calling Mosek 8.0.0.60: 28492 variables, 526 equality constraints

MOSEK Version 8.0.0.60 (Build date: 2017-3-1 13:09:33)
Copyright © MOSEK ApS, Denmark. WWW: mosek.com
Platform: Windows/64-X86

Problem
Name :
Objective sense : min
Type : CONIC (conic optimization problem)
Constraints : 526
Cones : 0
Scalar variables : 526
Matrix variables : 1
Integer variables : 0

Optimizer started.
Conic interior-point optimizer started.
Presolve started.
Linear dependency checker started.
Linear dependency checker terminated.
Eliminator - tries : 0 time : 0.00
Lin. dep. - tries : 1 time : 0.00
Lin. dep. - number : 0
Presolve terminated. Time: 0.02
Optimizer - threads : 4
Optimizer - solved problem : the primal
Optimizer - Constraints : 526
Optimizer - Cones : 1
Optimizer - Scalar variables : 527 conic : 55
Optimizer - Semi-definite variables: 1 scalarized : 27966
Factor - setup time : 0.03 dense det. time : 0.00
Factor - ML order time : 0.02 GP order time : 0.00
Factor - nonzeros before factor : 1.39e+005 after factor : 1.39e+005
Factor - dense dim. : 0 flops : 6.37e+007
ITE PFEAS DFEAS GFEAS PRSTATUS POBJ DOBJ MU TIME
0 2.0e+000 1.0e+000 1.0e+000 0.00e+000 0.000000000e+000 0.000000000e+000 1.0e+000 0.08
1 4.1e-001 2.1e-001 9.5e-002 -5.62e-001 -5.367173577e+001 -5.007263546e+001 2.1e-001 0.17
2 9.3e-002 4.7e-002 1.0e-002 -9.82e-001 -2.864147999e+002 -2.670895713e+002 4.7e-002 0.22
3 4.3e-002 2.2e-002 3.3e-003 -9.80e-001 -6.164326353e+002 -5.743814342e+002 2.2e-002 0.25
4 1.6e-002 7.8e-003 7.4e-004 -9.67e-001 -1.622493115e+003 -1.511709788e+003 7.8e-003 0.28
5 1.0e-002 5.3e-003 4.2e-004 -9.13e-001 -2.263922875e+003 -2.109595456e+003 5.3e-003 0.33
6 3.3e-003 1.7e-003 8.7e-005 -8.60e-001 -5.612962467e+003 -5.248070426e+003 1.7e-003 0.36
7 2.0e-003 1.0e-003 5.2e-005 -5.13e-001 -6.811677164e+003 -6.419791464e+003 1.0e-003 0.39
8 5.7e-004 2.9e-004 1.6e-005 -1.87e-001 -8.679107862e+003 -8.374770424e+003 2.9e-004 0.44
9 2.3e-004 1.1e-004 1.1e-005 6.50e-001 -6.050865842e+003 -5.938835721e+003 1.1e-004 0.47
10 5.4e-005 2.7e-005 6.1e-006 1.13e+000 -1.858101797e+003 -1.838579738e+003 2.7e-005 0.50
11 1.6e-005 8.2e-006 3.1e-006 9.83e-001 8.348984041e+001 9.055577170e+001 8.2e-006 0.53
12 6.5e-006 3.3e-006 2.0e-006 1.04e+000 9.140234458e+002 9.167672872e+002 3.3e-006 0.58
13 2.6e-006 1.3e-006 1.3e-006 1.09e+000 1.235220677e+003 1.236245382e+003 1.3e-006 0.61
14 9.3e-007 4.7e-007 7.6e-007 1.06e+000 1.367086387e+003 1.367462471e+003 4.7e-007 0.64
15 3.7e-007 1.9e-007 4.8e-007 1.03e+000 1.410511547e+003 1.410662579e+003 1.9e-007 0.69
16 1.5e-007 7.4e-008 3.0e-007 1.02e+000 1.428198945e+003 1.428259303e+003 7.4e-008 0.72
17 7.2e-008 3.6e-008 2.1e-007 1.01e+000 1.433941874e+003 1.433971258e+003 3.6e-008 0.75
18 2.5e-008 1.2e-008 1.2e-007 1.00e+000 1.437786263e+003 1.437796370e+003 1.2e-008 0.80
19 8.0e-009 4.0e-009 7.0e-008 1.00e+000 1.439040025e+003 1.439043332e+003 4.0e-009 0.83
20 1.7e-009 8.8e-010 3.2e-008 1.00e+000 1.439552273e+003 1.439552979e+003 8.6e-010 0.86
21 4.0e-010 7.3e-010 1.6e-008 1.00e+000 1.439654728e+003 1.439654891e+003 2.0e-010 0.91
22 8.8e-011 4.3e-010 7.3e-009 1.00e+000 1.439679208e+003 1.439679245e+003 4.4e-011 0.94
23 1.6e-011 8.2e-009 3.2e-009 1.00e+000 1.439684810e+003 1.439684817e+003 8.2e-012 0.98
Interior-point optimizer terminated. Time: 0.98.

Optimizer terminated. Time: 1.02

Interior-point solution summary
Problem status : PRIMAL_AND_DUAL_FEASIBLE
Solution status : OPTIMAL
Primal. obj: 1.4396848105e+003 nrm: 1e+003 Viol. con: 6e-004 var: 1e-004 barvar: 0e+000
Dual. obj: 1.4396848172e+003 nrm: 6e+003 Viol. con: 0e+000 var: 8e-009 barvar: 1e-011
Optimizer summary
Optimizer - time: 1.02
Interior-point - iterations : 23 time: 0.98
Basis identification - time: 0.00
Primal - iterations : 0 time: 0.00
Dual - iterations : 0 time: 0.00
Clean primal - iterations : 0 time: 0.00
Clean dual - iterations : 0 time: 0.00
Simplex - time: 0.00
Primal simplex - iterations : 0 time: 0.00
Dual simplex - iterations : 0 time: 0.00
Mixed integer - relaxations: 0 time: 0.00

Status: Solved
Optimal value (cvx_optval): +1439.68

computation time is around 117 seconds.

cvx .mtimes function take around 105 seconds. Here is the profile viewer information for this function:

image
image651×849 46.4 KB

AS you can see, around 99 seconds is about the self-time.

6

Thank you for your reply. Based on the code in the previous message, this is SDPT3 log:

Calling SDPT3 4.0: 28492 variables, 526 equality constraints

num. of constraints = 526
dim. of sdp var = 236, num. of sdp blk = 1
dim. of linear var = 472
dim. of free var = 54 *** convert ublk to lblk

SDPT3: Infeasible path-following algorithms

version predcorr gam expon scale_data
HKM 1 0.000 1 0
it pstep dstep pinfeas dinfeas gap prim-obj dual-obj cputime

0|0.000|0.000|2.0e+04|8.1e+04|5.6e+09| 4.191889e-05 0.000000e+00| 0:0:01| spchol 1 1
1|0.641|0.863|7.1e+03|1.1e+04|2.1e+09|-1.837380e+03 -8.524945e+06| 0:0:03| spchol 1 1
2|0.758|0.569|1.7e+03|4.8e+03|5.8e+08| 1.101166e+03 -1.039118e+07| 0:0:03| chol 1 1
3|0.232|0.728|1.3e+03|1.3e+03|5.8e+08| 1.860225e+03 -1.607828e+07| 0:0:03| chol 1 1
4|0.549|0.046|5.9e+02|1.3e+03|3.0e+08| 1.913134e+03 -1.644765e+07| 0:0:04| chol 1 1
5|0.553|0.360|2.6e+02|8.0e+02|1.7e+08| 2.401897e+03 -1.721941e+07| 0:0:04| chol 1 1
6|0.615|0.285|1.0e+02|5.7e+02|8.5e+07| 2.506307e+03 -1.647338e+07| 0:0:04| chol 1 1
7|0.672|0.332|3.3e+01|3.8e+02|4.1e+07| 2.754072e+03 -1.405402e+07| 0:0:04| chol 1 1
8|0.724|0.400|9.2e+00|2.3e+02|1.9e+07| 3.264459e+03 -1.012972e+07| 0:0:04| chol 1 1
9|0.768|0.421|2.1e+00|1.3e+02|9.5e+06| 4.216104e+03 -6.671176e+06| 0:0:04| chol 1 1
10|0.964|0.738|7.7e-02|3.5e+01|2.3e+06| 5.166830e+03 -2.099969e+06| 0:0:04| chol 1 1
11|0.986|0.983|1.1e-03|5.9e-01|4.1e+04| 5.223191e+03 -3.332341e+04| 0:0:04| chol 2 3
12|0.691|0.178|3.5e-04|4.9e-01|3.4e+04| 4.691425e+03 -2.726339e+04| 0:0:05| chol 2 3
13|0.230|0.256|2.7e-04|3.6e-01|2.7e+04| 4.514930e+03 -2.058853e+04| 0:0:05| chol 2 3
14|0.504|0.105|1.3e-04|3.2e-01|2.4e+04| 4.053380e+03 -1.855431e+04| 0:0:05| chol 2 4
15|0.488|0.328|6.7e-05|2.2e-01|1.8e+04| 3.613876e+03 -1.274078e+04| 0:0:05| chol 2 4
16|0.603|0.219|2.7e-05|1.7e-01|1.4e+04| 2.976939e+03 -9.965290e+03| 0:0:05| chol 2 3
17|0.622|0.295|1.0e-05|1.2e-01|9.9e+03| 2.407845e+03 -6.843247e+03| 0:0:05| chol 3 4
18|0.707|0.380|3.0e-06|7.4e-02|6.0e+03| 1.903096e+03 -3.810983e+03| 0:0:05| chol 3 5
19|0.497|0.188|1.5e-06|6.0e-02|4.9e+03| 1.811203e+03 -2.880644e+03| 0:0:05| chol 3 5
20|0.650|0.463|5.2e-07|3.2e-02|2.7e+03| 1.663721e+03 -9.406696e+02| 0:0:06| chol 4 5
21|0.879|0.461|6.4e-08|1.7e-02|1.4e+03| 1.512344e+03 1.479144e+02| 0:0:06| chol 4 7
22|1.000|0.131|1.6e-08|1.5e-02|1.3e+03| 1.517184e+03 3.102598e+02| 0:0:06| chol 4 5
23|1.000|0.597|7.5e-09|6.1e-03|5.4e+02| 1.482056e+03 9.707428e+02| 0:0:06| chol 4 6
24|1.000|0.231|7.1e-09|4.7e-03|4.3e+02| 1.474117e+03 1.073957e+03| 0:0:06| chol 4 6
25|1.000|0.264|6.7e-09|3.5e-03|3.4e+02| 1.471262e+03 1.164511e+03| 0:0:06| chol 4 5
26|1.000|0.374|3.9e-09|2.2e-03|2.2e+02| 1.458634e+03 1.262378e+03| 0:0:06| chol 4 3
27|1.000|0.324|2.7e-09|1.5e-03|1.5e+02| 1.452477e+03 1.316947e+03| 0:0:06| chol 4 5
28|1.000|0.349|1.7e-09|9.5e-04|9.8e+01| 1.447727e+03 1.357773e+03| 0:0:07| chol 3 3
29|1.000|0.365|9.9e-10|6.1e-04|6.3e+01| 1.444360e+03 1.386485e+03| 0:0:07| chol 3 3
30|1.000|0.362|9.4e-10|3.9e-04|4.0e+01| 1.442400e+03 1.405088e+03| 0:0:07| chol 4 3
31|1.000|0.363|3.3e-10|2.5e-04|2.6e+01| 1.441252e+03 1.417256e+03| 0:0:07| chol 3 3
32|1.000|0.374|1.9e-10|1.5e-04|1.6e+01| 1.440548e+03 1.425426e+03| 0:0:07| chol 3 3
33|1.000|0.387|9.0e-11|9.4e-05|9.8e+00| 1.440131e+03 1.430836e+03| 0:0:07| chol 3 3
34|1.000|0.388|4.6e-11|5.8e-05|5.9e+00| 1.439913e+03 1.434208e+03| 0:0:07| chol 3 3
35|1.000|0.402|2.4e-11|3.5e-05|3.5e+00| 1.439795e+03 1.436383e+03| 0:0:07| chol 3 3
36|1.000|0.415|1.0e-11|2.0e-05|2.0e+00| 1.439735e+03 1.437741e+03| 0:0:07| chol 3 3
37|1.000|0.377|6.2e-12|1.3e-05|1.3e+00| 1.439712e+03 1.438468e+03| 0:0:08| chol 3 3
38|1.000|0.419|3.5e-12|7.3e-06|7.4e-01| 1.439698e+03 1.438975e+03| 0:0:08| chol 3 3
39|1.000|0.410|3.8e-12|4.3e-06|4.3e-01| 1.439691e+03 1.439265e+03| 0:0:08| chol 3 3
40|1.000|0.438|2.8e-12|2.4e-06|2.4e-01| 1.439688e+03 1.439449e+03| 0:0:08| chol 4 3
41|1.000|0.492|3.2e-12|1.2e-06|1.2e-01| 1.439687e+03 1.439565e+03| 0:0:08| chol 3 3
42|1.000|0.541|3.6e-12|5.7e-07|5.6e-02| 1.439686e+03 1.439631e+03| 0:0:08| chol 4 4
43|1.000|0.612|2.8e-12|3.0e-06|2.5e-02| 1.439686e+03 1.439665e+03| 0:0:08| chol 8 8
44|1.000|0.974|5.2e-12|6.6e-06|2.3e-03| 1.439686e+03 1.439686e+03| 0:0:08| chol 10 10
45|0.932|0.973|6.2e-11|9.8e-07|1.1e-04| 1.439686e+03 1.439686e+03| 0:0:09| chol
linsysolve: Schur complement matrix not positive definite
switch to LU factor. lu 11 ^ 1
46|1.000|0.888|4.1e-01|4.8e-08|9.0e-04| 1.573593e+03 1.439686e+03| 0:0:09| lu 11 ^ 3
47|0.968|0.748|9.4e-01|3.8e-07|5.7e-05| 1.594758e+03 1.439686e+03| 0:0:09| lu 5 1
48|0.126|0.699|8.2e-01|8.7e-08|4.8e-04| 1.575142e+03 1.439683e+03| 0:0:09| lu 3 1
49|0.960|0.302|3.3e-02|1.2e-07|2.2e-02| 1.445134e+03 1.439664e+03| 0:0:09| lu 5 1
50|1.000|0.893|8.6e-12|9.4e-06|3.7e-03| 1.439686e+03 1.439684e+03| 0:0:09| lu 5 1
51|0.969|0.984|8.6e-12|1.6e-06|1.6e-04| 1.439686e+03 1.439686e+03| 0:0:10| lu 11 ^ 5
52|0.922|0.891|3.0e+00|6.9e-08|1.5e-04| 1.538799e+03 1.439686e+03| 0:0:10| lu 15 1
53|0.620|0.844|1.1e+00|1.0e-07|1.0e-03| 1.477363e+03 1.439685e+03| 0:0:10| lu 30 30
54|0.946|0.734|6.2e-02|4.5e-07|6.7e-04| 1.441735e+03 1.439686e+03| 0:0:10| lu 11 30
55|0.000|0.000|6.6e-02|3.4e-07|6.5e-04| 1.441778e+03 1.439686e+03| 0:0:10| lu 30 30
56|0.978|0.932|1.4e+00|2.8e-07|9.2e-05| 1.119489e+03 1.439686e+03| 0:0:10| lu 30 30
57|0.944|0.438|1.1e+01|4.3e-08|1.7e-05| 7.847928e+02 1.439686e+03| 0:0:11| lu 30 30
58|0.005|0.040|4.0e+01|1.3e-08|1.5e-05| 7.821195e+02 1.439686e+03| 0:0:11| lu 30 30
59|0.000|0.005|5.8e+01|1.2e-08|1.5e-05| 7.825654e+02 1.439686e+03| 0:0:11| lu 30 30
60|0.000|0.002|6.7e+01|1.2e-08|1.5e-05| 7.825917e+02 1.439686e+03| 0:0:11| lu 30 30
61|0.000|0.000|7.6e+01|1.2e-08|1.5e-05| 7.825941e+02 1.439686e+03| 0:0:11| lu 30 30
stop: Z not positive definite
62|0.000|0.000|7.6e+01|1.2e-08|1.5e-05| 7.825941e+02 1.439686e+03| 0:0:12|

number of iterations = 62
primal objective value = 1.43968610e+03
dual objective value = 1.43968609e+03
gap := trace(XZ) = 1.11e-04
relative gap = 3.85e-08
actual relative gap = 5.74e-09
rel. primal infeas (scaled problem) = 6.24e-11
rel. dual " " " = 9.76e-07
rel. primal infeas (unscaled problem) = 0.00e+00
rel. dual " " " = 0.00e+00
norm(X), norm(y), norm(Z) = 3.6e+03, 2.5e+06, 2.5e+06
norm(A), norm(b), norm© = 2.2e+05, 3.6e+03, 1.1e+01
Total CPU time (secs) = 11.57
CPU time per iteration = 0.19
termination code = -3
DIMACS: 2.2e-10 0.0e+00 5.6e-06 0.0e+00 5.7e-09 3.8e-08

Status: Inaccurate/Solved
Optimal value (cvx_optval): +1439.69

computation time is around 102 seconds which is faster than mosek. However, around 87 seconds are spent on the cvx.mtimes line 203.

7

It looks like Mosek only took a second to solve the problem. So if you are going to get much improvement in end to end wall clock time, it will have to be by reducing CVX model processing time. Can you vectorize constraint generation rather than using for loops?

8

Y_k and Y_k_bar in my code come from a handle function and also there is a trace involved in the constraints. I’m not sure whether I can use vectorize constraint instead of for loop. I need to check it.

9

Update: vectorizing the constraints does not have any positive impacts on the computation burden.

10

Can you show us the vectorized code?

12 
for i=1:length(Nbus)
    Y_kmega{i}=Y_k(i);
end
 Y_kmega=reshape(Y_kmega,[length(Nbus),1]);

 Y_kmega=cell2mat(Y_kmega);

  for obj_gen_idx=1:length(Gloc)
      obj_gen1=Gloc(obj_gen_idx);
      Y_kmega1{obj_gen_idx}=Y_k(obj_gen1);
end

 Y_kmega1=reshape(Y_kmega1,[length(Gloc),1]);

Y_kmega1=cell2mat(Y_kmega1);


 cvx_begin sdp
    cvx_solver mosek

    variable W(2*length(Nbus),2*length(Nbus)) 
    variable obj_gen
  
  minimize obj_gen
  subject  to

%% constraint related to upper and lower bound for generated active ad reactive power

     obj_gen==100*trace(Y_kmega1*repmat(W,1,length(Gloc)));
    100*sum(reshape(diag((Y_kmega*repmat(W,1,length(Nbus)))),2*Nbus,[]))'<=Pmax-PD;
    100*sum(reshape(diag((Y_kmega*repmat(W,1,length(Nbus)))),2*Nbus,[]))'>=Pmin-PD;

% 100*trace(Y_k_bar(g_idx)W)<=Qmax(g_idx)-QD(g_idx);
% 100trace(Y_k_bar(g_idx)*W)>=Qmin(g_idx)-QD(g_idx);

%% semidefinite Constraint

      W>=0;

    cvx_end

Unfortunately, it took a lot of time to solve that I canceled it (a lot more than for loop code-out of memory). I don’t know if there is a better way to vectorized the constraints.