Using small numbers (~1e-20) in CVX

[nikku]

I am currently working on a SDP problem as follows:

%% define problem
% init
cvx_begin sdp quiet
cvx_precision high

% define variables
variables r12 r13 r23

% intermediate variables
r33 = s12 - r12 + r13 + r23;
r22 = s22 - r33 + 2*r23;
r11 = s11 - r33 + 2*r13;

R = [r11,r12,r13;
     r12,r22,r23;
     r13,r23,r33];
 
F = (r12^2 + r13^2 + r23^2);

% objective
minimize(F)

% constraints
subject to:
R >= 0


cvx_end

where s_{11},s_{12},s_{22} are provided. The problem is solved fine when the values for s_{11},s_{12},s_{22} are on a scale of 1e-10 or more. But, when I use the desired values (on a scale of 1e-20), the solution is strange (provides a feasible solution that isn’t minimized very much).

I think think this is a numerical issue regarding precision. So my question is can this problem be scaled to obtain a cleaner solution?

Below are the results for three sets of s_{11},s_{12},s_{22} values to illustrate what I mean

SET 1

% external parameters
alph = 1e20;
s11 = 1.6227e-20 * alph;
s12 = 2.3962e-22 * alph;
s22 = 1.3187e-20 * alph;

results in

Calling SDPT3 4.0: 15 variables, 6 equality constraints
   For improved efficiency, SDPT3 is solving the dual problem.
------------------------------------------------------------

 num. of constraints =  6
 dim. of sdp    var  =  9,   num. of sdp  blk  =  4
*******************************************************************
   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|8.5e+00|7.6e+00|9.0e+02| 5.917438e+01  0.000000e+00| 0:0:00| chol  1  1 
 1|1.000|1.000|3.0e-06|8.1e-02|7.1e+01| 4.179353e+01 -2.498514e+01| 0:0:00| chol  1  1 
 2|0.950|0.943|8.2e-07|1.2e-02|4.0e+00| 2.754911e+00 -1.125021e+00| 0:0:00| chol  1  1 
 3|1.000|1.000|1.3e-07|8.1e-04|1.2e+00| 7.398372e-01 -4.271277e-01| 0:0:00| chol  1  1 
 4|0.922|0.917|2.1e-08|1.4e-04|1.1e-01| 6.988350e-02 -4.350303e-02| 0:0:00| chol  1  1 
 5|1.000|1.000|2.5e-09|8.1e-06|4.2e-02| 2.702564e-02 -1.487091e-02| 0:0:00| chol  1  1 
 6|0.938|0.938|5.2e-10|1.3e-06|2.7e-03| 1.981539e-03 -7.339812e-04| 0:0:00| chol  1  1 
 7|1.000|1.000|2.0e-09|8.1e-08|7.0e-04| 4.927488e-04 -2.069753e-04| 0:0:00| chol  1  1 
 8|0.946|0.945|4.9e-10|1.2e-08|4.5e-05| 3.861217e-05 -6.345292e-06| 0:0:00| chol  1  1 
 9|0.983|0.991|3.4e-10|2.1e-10|1.2e-06| 1.047654e-06 -1.594209e-07| 0:0:00| chol  1  1 
10|0.992|1.000|2.7e-12|6.8e-11|3.0e-08| 2.150749e-08 -7.745064e-09| 0:0:00| chol  1  1 
11|1.000|1.000|5.7e-17|1.0e-12|1.1e-09| 6.681392e-10 -4.184860e-10| 0:0:00| chol  1  1 
12|1.000|1.000|9.1e-17|1.0e-12|4.1e-11| 2.359865e-11 -1.381214e-11| 0:0:00| chol  1  1 
13|1.000|1.000|7.2e-23|1.8e-14|7.2e-13| 4.117973e-13 -2.402632e-13| 0:0:00|
  stop: max(relative gap, infeasibilities) < 1.82e-12
-------------------------------------------------------------------
 number of iterations   = 13
 primal objective value =  4.11797339e-13
 dual   objective value = -2.40263194e-13
 gap := trace(XZ)       = 7.17e-13
 relative gap           = 7.17e-13
 actual relative gap    = 6.52e-13
 rel. primal infeas (scaled problem)   = 7.19e-23
 rel. dual     "        "       "      = 1.77e-14
 rel. primal infeas (unscaled problem) = 0.00e+00
 rel. dual     "        "       "      = 0.00e+00
 norm(X), norm(y), norm(Z) = 1.7e+00, 5.6e-08, 2.7e+00
 norm(A), norm(b), norm(C) = 5.9e+00, 2.7e+00, 3.7e+00
 Total CPU time (secs)  = 0.11  
 CPU time per iteration = 0.01  
 termination code       =  0
 DIMACS: 9.8e-23  0.0e+00  2.5e-14  0.0e+00  6.5e-13  7.2e-13
-------------------------------------------------------------------
 
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +2.40263e-13

R =

    1.5987    0.0000   -0.0000
    0.0000    1.2947   -0.0000
   -0.0000   -0.0000    0.0240

Eigenvalues of R =  0.0240    1.2947    1.5987

SET 2

% external parameters
alph = 1e15;
s11 = 1.6227e-20 * alph;
s12 = 2.3962e-22 * alph;
s22 = 1.3187e-20 * alph;

results in

Calling SDPT3 4.0: 15 variables, 6 equality constraints
   For improved efficiency, SDPT3 is solving the dual problem.
------------------------------------------------------------

 num. of constraints =  6
 dim. of sdp    var  =  9,   num. of sdp  blk  =  4
*******************************************************************
   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|8.5e+00|1.1e+01|9.0e+02| 3.000029e+01  0.000000e+00| 0:0:00| chol  1  1 
 1|1.000|0.999|3.2e-06|1.1e-01|5.5e+01| 2.499765e+01 -2.518871e+01| 0:0:00| chol  1  1 
 2|0.949|0.942|4.9e-06|1.6e-02|3.3e+00| 1.551908e+00 -1.423246e+00| 0:0:00| chol  1  1 
 3|1.000|1.000|1.9e-06|1.0e-03|7.4e-01| 3.309290e-01 -3.833775e-01| 0:0:00| chol  1  1 
 4|0.922|0.908|1.6e-06|1.8e-04|8.6e-02| 4.136656e-02 -3.964401e-02| 0:0:00| chol  1  1 
 5|1.000|1.000|6.6e-08|1.0e-05|3.1e-02| 1.527952e-02 -1.516405e-02| 0:0:00| chol  1  1 
 6|0.923|0.924|4.1e-08|1.7e-06|2.8e-03| 1.457818e-03 -1.252489e-03| 0:0:00| chol  1  1 
 7|1.000|1.000|7.5e-09|1.1e-07|1.0e-03| 6.064261e-04 -4.010489e-04| 0:0:00| chol  1  1 
 8|0.921|0.921|1.5e-09|1.9e-08|9.2e-05| 5.504201e-05 -3.664925e-05| 0:0:00| chol  1  1 
 9|1.000|1.000|3.0e-09|3.0e-10|3.7e-05| 2.246353e-05 -1.430350e-05| 0:0:00| chol  1  1 
10|0.928|0.928|2.2e-10|4.6e-10|2.9e-06| 1.759897e-06 -1.112182e-06| 0:0:00| chol  1  1 
11|1.000|1.000|1.1e-14|4.3e-11|1.1e-06| 6.873935e-07 -4.324245e-07| 0:0:00| chol  1  1 
12|1.000|1.000|2.0e-15|1.0e-12|1.3e-07| 7.875050e-08 -4.882515e-08| 0:0:00| chol  1  1 
13|1.000|1.000|2.2e-16|1.0e-12|2.2e-08| 1.362100e-08 -8.134173e-09| 0:0:00| chol  1  1 
14|1.000|1.000|2.9e-16|1.0e-12|3.4e-09| 2.207793e-09 -1.227905e-09| 0:0:00| chol  1  1 
15|1.000|1.000|1.1e-16|1.0e-12|5.6e-10| 3.725307e-10 -1.884100e-10| 0:0:00| chol  1  1 
16|1.000|1.000|1.2e-16|1.0e-12|9.8e-11| 6.155822e-11 -3.392790e-11| 0:0:00| chol  1  1 
17|1.000|1.000|2.1e-16|6.0e-13|1.8e-11| 1.036226e-11 -6.028238e-12| 0:0:00| chol  1  1 
18|1.000|1.000|2.6e-16|7.9e-14|2.4e-12| 1.489879e-12 -6.591993e-13| 0:0:00| chol  1  1 
19|1.000|1.000|1.6e-16|1.1e-14|3.3e-13| 2.211559e-13 -8.032774e-14| 0:0:00|
  stop: max(relative gap, infeasibilities) < 1.82e-12
-------------------------------------------------------------------
 number of iterations   = 19
 primal objective value =  2.21155907e-13
 dual   objective value = -8.03277438e-14
 gap := trace(XZ)       = 3.32e-13
 relative gap           = 3.32e-13
 actual relative gap    = 3.01e-13
 rel. primal infeas (scaled problem)   = 1.63e-16
 rel. dual     "        "       "      = 1.10e-14
 rel. primal infeas (unscaled problem) = 0.00e+00
 rel. dual     "        "       "      = 0.00e+00
 norm(X), norm(y), norm(Z) = 1.7e+00, 1.6e-07, 1.7e+00
 norm(A), norm(b), norm(C) = 5.9e+00, 2.7e+00, 2.7e+00
 Total CPU time (secs)  = 0.13  
 CPU time per iteration = 0.01  
 termination code       =  0
 DIMACS: 2.2e-16  0.0e+00  1.5e-14  0.0e+00  3.0e-13  3.3e-13
-------------------------------------------------------------------
 
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +8.03277e-14

R =
   1.0e-04 *

    0.1590   -0.0009    0.0009
   -0.0009    0.1286    0.0009
    0.0009    0.0009    0.0052

Eigenvalues of R =  1.0e-04 *
    0.0051    0.1286    0.1590

SET 3

% external parameters
alph = 1e10;
s11 = 1.6227e-20 * alph;
s12 = 2.3962e-22 * alph;
s22 = 1.3187e-20 * alph;

results in

Calling SDPT3 4.0: 15 variables, 6 equality constraints
   For improved efficiency, SDPT3 is solving the dual problem.
------------------------------------------------------------

 num. of constraints =  6
 dim. of sdp    var  =  9,   num. of sdp  blk  =  4
*******************************************************************
   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|8.5e+00|1.1e+01|9.0e+02| 3.000000e+01  0.000000e+00| 0:0:00| chol  1  1 
 1|1.000|0.999|3.2e-06|1.1e-01|5.5e+01| 2.499743e+01 -2.518870e+01| 0:0:00| chol  1  1 
 2|0.949|0.942|4.9e-06|1.6e-02|3.3e+00| 1.551742e+00 -1.423247e+00| 0:0:00| chol  1  1 
 3|1.000|1.000|1.9e-06|1.0e-03|7.4e-01| 3.306158e-01 -3.833480e-01| 0:0:00| chol  1  1 
 4|0.922|0.908|1.7e-06|1.8e-04|8.6e-02| 4.108610e-02 -3.965060e-02| 0:0:00| chol  1  1 
 5|1.000|1.000|6.9e-08|1.0e-05|3.0e-02| 1.454036e-02 -1.510419e-02| 0:0:00| chol  1  1 
 6|0.922|0.923|7.1e-08|1.7e-06|2.6e-03| 1.256959e-03 -1.254484e-03| 0:0:00| chol  1  1 
 7|1.000|1.000|2.5e-08|1.1e-07|7.9e-04| 3.896620e-04 -3.817790e-04| 0:0:00| chol  1  1 
 8|0.916|0.916|1.5e-08|2.4e-08|7.8e-05| 3.673919e-05 -3.670699e-05| 0:0:00| chol  1  1 
 9|1.000|0.954|3.0e-09|4.2e-09|2.8e-05| 1.269306e-05 -1.374350e-05| 0:0:00| chol  1  1 
10|0.918|0.923|2.7e-09|8.6e-10|2.8e-06| 1.225321e-06 -1.203870e-06| 0:0:00| chol  1  1 
11|1.000|1.000|5.1e-09|1.5e-10|1.1e-06| 5.366496e-07 -4.304362e-07| 0:0:00| chol  1  1 
12|1.000|1.000|1.4e-08|3.8e-11|1.7e-07| 9.332970e-08 -6.165906e-08| 0:0:00| chol  1  1 
13|1.000|1.000|8.8e-09|3.0e-11|3.7e-08| 2.129104e-08 -1.244525e-08| 0:0:00| chol  1  1 
14|1.000|1.000|1.1e-09|2.4e-11|6.9e-09| 3.980698e-09 -2.329622e-09| 0:0:00| chol  1  1 
15|1.000|1.000|4.8e-11|1.6e-11|1.2e-09| 7.043804e-10 -4.010195e-10| 0:0:00| chol  1  1 
16|1.000|1.000|8.0e-12|5.0e-12|1.8e-10| 1.056404e-10 -5.954175e-11| 0:0:00| chol  1  1 
17|1.000|1.000|9.1e-13|7.3e-13|2.2e-11| 1.326695e-11 -7.092128e-12| 0:0:00| chol  1  1 
18|1.000|1.000|8.1e-14|1.1e-13|3.2e-12| 1.914561e-12 -1.017984e-12| 0:0:00| chol  1  1 
19|1.000|1.000|1.1e-13|1.7e-14|5.2e-13| 3.084508e-13 -1.641854e-13| 0:0:00|
  stop: max(relative gap, infeasibilities) < 1.82e-12
-------------------------------------------------------------------
 number of iterations   = 19
 primal objective value =  3.08450828e-13
 dual   objective value = -1.64185379e-13
 gap := trace(XZ)       = 5.20e-13
 relative gap           = 5.20e-13
 actual relative gap    = 4.73e-13
 rel. primal infeas (scaled problem)   = 1.08e-13
 rel. dual     "        "       "      = 1.73e-14
 rel. primal infeas (unscaled problem) = 0.00e+00
 rel. dual     "        "       "      = 0.00e+00
 norm(X), norm(y), norm(Z) = 1.7e+00, 2.8e-07, 1.7e+00
 norm(A), norm(b), norm(C) = 5.9e+00, 2.7e+00, 2.7e+00
 Total CPU time (secs)  = 0.10  
 CPU time per iteration = 0.01  
 termination code       =  0
 DIMACS: 1.5e-13  0.0e+00  2.4e-14  0.0e+00  4.7e-13  5.2e-13
-------------------------------------------------------------------
 
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +1.64185e-13

R =
   1.0e-06 *

    0.1631    0.1629    0.1630
    0.1629    0.1631    0.1630
    0.1630    0.1630    0.1630

Eigenvalues of R =     1.0e-06 *
    0.0000    0.0001    0.4890

Notice that the first two sets have solutions for R where r_{12}, r_{13}, r_{23} are clearly minimized. You can see it by the fact that they are smaller than the diagonal terms.

But for the third set, r_{12}, r_{13}, r_{23} are close to the diagonal terms and seem like they were not minimized very much.

[Mark_L_Stone]

Double precision solvers, such as those called by CVX, don’t like numbers such as 1e-20. Use of 1e-10 is still rather dubious. Try to keep all non-zero input data within a small number of orders of magnitude of one.

Due to its greater numerical robustness, you might get away with more if Mosek is used, rather than SDPT3 or SeDuMi. But you should still try to avoid badly scaled input data. If using Mosrk, you should pay attention to warning messages about large or nearly zero input data.

[nikku]

@Mark_L_Stone thanks for the reply! Yeah, the input data scale is unreasonable, I agree. Can you point me to any resources for how to rescale the problem?

I find that there isn’t easily accessible information on the internet about how to rescale optimization problems. Maybe because scaling isn’t regarding theory but rather numerical techniques. So it seems like there are some tricks of the trade, which I am unaware of.

[Mark_L_Stone]

Try to choose units, for instance picowatts instead of watts, to improve scaling.

You can also read
sections 7.2 “Avoiding ill-posed problems” and 7.3 “Scaling” of 7 Practical optimization — MOSEK Modeling Cookbook 3.3.0

and

https://docs.mosek.com/latest/capi/debugging-numerical.html

[Michal_Adamaszek]

Here is something about rescaling I wrote in another thread.

How would you enter the following expression into cvx?f（x,y)=x^3/y^2 Functions like pow don't seem to work as a whole for cvx to recognize convexity

There are two energy drinks. One gives 10^{-18} calories per litre, the other 1 calorie per litre. You want to find some combination of them which maximizes the total number of calories 10^{-18}x+y, subject to some other constraints, where x,y denote how many litres of each drink you want to intake. There are two possibilities. Either you expect x,y to be of reasonable magnitude for a human being to drink, say somewhere about 1 litre. In that case x completely doesn’t matter, including it is w…