getdual in conic interface to Mosek throws BoundsError

[chriscoey]

Here is an example of the error I get when I call getdual:

ERROR: LoadError: BoundsError: attempt to access 2-element Array{Array{Float64,1},1}:
 [0.0666205723355703,0.04302634696803116,-0.25810960668167754,0.027788215544164746,-0.16669795651893662,1.0000000011590489]
 [0.027774525875638566,0.09856627300800122,-0.16665690730809762,0.3497921239111226,-0.5914322621022672,1.0000000021142892]
  at index [3]
 in getdual at /Users/chris/.julia/v0.4/Mosek/src/MosekConicInterface.jl:526

Unfortunately, I don't understand what's going on here at line 526:

Float64[if     m.conslack[i] == 0 y[i]
            elseif m.conslack[i] >  0 snx[m.conslack[i]]
            else                      bars[-m.conslack[i]][m.barconij[i]]
            end
            for i in 1:m.numcon ]

[mlubin]

Could you provide an example script that reproduces this error ?

[chriscoey]

Here is a script (I can produce more if it is helpful):

c = [0.0,0.0,1.0,0.0,0.0,0.0,0.0]
b = [0.0,10.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,0.0]

I = [1,12,1,21,1,2,3,7,8,10,13,16,17,19,22,2,4,7,16,2,5,7,16,2,6,10,19]
J = [1,1,2,2,3,4,4,4,4,4,4,4,4,4,4,5,5,5,5,6,6,6,6,7,7,7,7]
V = [-1.0,-1.0,-1.0,-1.0,1.0,1.0,-1.0,-0.44999999999999996,0.45,-0.4500000000000001,0.,-0.44999999999999996,0.45,-0.4500000000000001,0.,1.0,-1.0,-0.9000000000000001,-0.9000000000000001,1.0,-1.0,-0.22500000000000003,-0.22500000000000003,1.0,-1.0,-0.22500000000000003,-0.22500000000000003]
A = sparse(I, J, V, length(b), length(c))

cone_con = [(:Zero,1:1),(:NonNeg,2:2),(:NonNeg,3:6),(:SDP,7:12),(:Zero,13:15),(:SDP,16:21),(:Zero,22:24)]
cone_var = [(:Free,1:7)]

using MathProgBase, Mosek

m = MathProgBase.ConicModel(MosekSolver())
MathProgBase.loadproblem!(m, c, A, b, cone_con, cone_var)
MathProgBase.optimize!(m)

@show MathProgBase.status(m)
@show MathProgBase.getsolution(m)
@show MathProgBase.getdual(m)

Here is the output I get:

julia> MathProgBase.optimize!(m)
Computer
  Platform               : MACOSX/64-X86   

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

Optimizer started.
Conic interior-point optimizer started.
Presolve started.
Linear dependency checker started.
Linear dependency checker terminated.
Eliminator started.
Total number of eliminations : 1
Eliminator terminated.
Eliminator - tries                  : 1                 time                   : 0.00            
Eliminator - elim's                 : 1               
Lin. dep.  - tries                  : 1                 time                   : 0.00            
Lin. dep.  - number                 : 0               
Presolve terminated. Time: 0.00    
Optimizer  - threads                : 8               
Optimizer  - solved problem         : the primal      
Optimizer  - Constraints            : 13
Optimizer  - Cones                  : 0
Optimizer  - Scalar variables       : 8                 conic                  : 0               
Optimizer  - Semi-definite variables: 2                 scalarized             : 12              
Factor     - setup time             : 0.00              dense det. time        : 0.00            
Factor     - ML order time          : 0.00              GP order time          : 0.00            
Factor     - nonzeros before factor : 59                after factor           : 63              
Factor     - dense dim.             : 0                 flops                  : 1.01e+03        
ITE PFEAS    DFEAS    GFEAS    PRSTATUS   POBJ              DOBJ              MU       TIME  
0   5.0e+00  1.0e+00  1.0e+00  0.00e+00   0.000000000e+00   0.000000000e+00   1.0e+00  0.00  
1   1.3e+00  2.6e-01  2.6e-01  -1.45e-01  7.446299248e-01   1.133476135e+00   2.6e-01  0.00  
2   1.8e-01  3.6e-02  3.6e-02  5.79e-01   1.297888733e+00   1.369246549e+00   3.6e-02  0.00  
3   2.4e-02  4.7e-03  4.7e-03  8.91e-01   9.564969913e-01   9.704884077e-01   4.7e-03  0.00  
4   3.4e-03  6.8e-04  6.8e-04  8.77e-01   8.651595034e-01   8.677722067e-01   6.8e-04  0.00  
5   2.7e-04  5.4e-05  5.4e-05  9.69e-01   8.507793221e-01   8.509928746e-01   5.4e-05  0.00  
6   1.6e-05  3.3e-06  3.3e-06  9.97e-01   8.496036890e-01   8.496167452e-01   3.3e-06  0.00  
7   9.2e-07  1.8e-07  1.8e-07  1.00e+00   8.495321450e-01   8.495328766e-01   1.8e-07  0.00  
8   8.8e-09  1.8e-09  1.8e-09  1.00e+00   8.495279655e-01   8.495279725e-01   1.8e-09  0.00  
Interior-point optimizer terminated. Time: 0.00. 

Optimizer terminated. Time: 0.02    
0

julia> @show MathProgBase.status(m)
MathProgBase.status(m) = :Optimal
:Optimal

julia> @show MathProgBase.getsolution(m)
MathProgBase.getsolution(m) = [0.2581288811960397,0.591399084272713,0.8495279654687528,2.965654193424076,3.779037063878916,3.49217819844456e-7,3.255308073884936]
7-element Array{Float64,1}:
 0.258129  
 0.591399  
 0.849528  
 2.96565   
 3.77904   
 3.49218e-7
 3.25531   

julia> @show MathProgBase.getdual(m)
ERROR: BoundsError: attempt to access 2-element Array{Array{Float64,1},1}:
 [0.06662544466703094,0.04302513616541794,-0.25811904947274705,0.027784617540638323,-0.16668717064440128,0.9999999997848541]
 [0.027780822124551928,0.09858424790224177,-0.1666757940939465,0.3498404168267607,-0.5914730859020215,0.9999999999546839]   
  at index [3]
 in getdual at /Users/chris/.julia/v0.4/Mosek/src/MosekConicInterface.jl:526

[chriscoey]

for comparison, here is the output I get from SCS solver:

julia> MathProgBase.optimize!(m)
----------------------------------------------------------------------------
    SCS v1.1.8 - Splitting Conic Solver
    (c) Brendan O'Donoghue, Stanford University, 2012-2015
----------------------------------------------------------------------------
Lin-sys: sparse-direct, nnz in A = 25
eps = 1.00e-04, alpha = 1.80, max_iters = 20000, normalize = 1, scale = 5.00
Variables n = 7, constraints m = 24
Cones:  primal zero / dual free vars: 7
    linear vars: 5
    sd vars: 12, sd blks: 2
Setup time: 1.91e-03s
----------------------------------------------------------------------------
 Iter | pri res | dua res | rel gap | pri obj | dua obj | kap/tau | time (s)
----------------------------------------------------------------------------
     0|      inf       inf       nan      -inf       inf       inf  5.12e-04 
   100| 4.04e-04  2.37e-02  3.95e-03  9.03e-01  9.14e-01  0.00e+00  1.16e-03 
   200| 3.17e-04  1.62e-02  4.02e-03  8.83e-01  8.94e-01  0.00e+00  1.79e-03 
   300| 2.56e-04  1.19e-02  3.48e-03  8.72e-01  8.81e-01  0.00e+00  2.41e-03 
   400| 2.10e-04  9.02e-03  3.00e-03  8.64e-01  8.73e-01  0.00e+00  3.03e-03 
   500| 1.74e-04  7.04e-03  2.57e-03  8.60e-01  8.67e-01  0.00e+00  3.65e-03 
   600| 1.45e-04  5.59e-03  2.19e-03  8.57e-01  8.63e-01  0.00e+00  4.27e-03 
   700| 1.21e-04  4.49e-03  1.86e-03  8.54e-01  8.59e-01  0.00e+00  4.89e-03 
   800| 1.01e-04  3.64e-03  1.58e-03  8.53e-01  8.57e-01  0.00e+00  5.51e-03 
   900| 8.46e-05  2.97e-03  1.34e-03  8.52e-01  8.56e-01  0.00e+00  6.12e-03 
  1000| 7.09e-05  2.44e-03  1.13e-03  8.51e-01  8.54e-01  0.00e+00  6.74e-03 
  1100| 5.94e-05  2.01e-03  9.58e-04  8.51e-01  8.53e-01  0.00e+00  7.36e-03 
  1200| 4.98e-05  1.66e-03  8.07e-04  8.50e-01  8.53e-01  0.00e+00  7.98e-03 
  1300| 4.17e-05  1.37e-03  6.80e-04  8.50e-01  8.52e-01  0.00e+00  8.59e-03 
  1400| 3.50e-05  1.14e-03  5.72e-04  8.50e-01  8.52e-01  0.00e+00  9.21e-03 
  1500| 2.93e-05  9.47e-04  4.80e-04  8.50e-01  8.51e-01  0.00e+00  9.82e-03 
  1600| 2.45e-05  7.88e-04  4.03e-04  8.50e-01  8.51e-01  0.00e+00  1.05e-02 
  1700| 2.05e-05  6.56e-04  3.39e-04  8.50e-01  8.51e-01  0.00e+00  1.11e-02 
  1800| 1.72e-05  5.47e-04  2.84e-04  8.50e-01  8.50e-01  0.00e+00  1.17e-02 
  1900| 1.44e-05  4.56e-04  2.38e-04  8.50e-01  8.50e-01  0.00e+00  1.23e-02 
  2000| 1.20e-05  3.80e-04  2.00e-04  8.50e-01  8.50e-01  0.00e+00  1.29e-02 
  2100| 1.01e-05  3.17e-04  1.67e-04  8.50e-01  8.50e-01  0.00e+00  1.35e-02 
  2200| 8.44e-06  2.65e-04  1.40e-04  8.50e-01  8.50e-01  0.00e+00  1.41e-02 
  2300| 7.06e-06  2.21e-04  1.17e-04  8.50e-01  8.50e-01  0.00e+00  1.48e-02 
  2400| 5.91e-06  1.85e-04  9.81e-05  8.50e-01  8.50e-01  0.00e+00  1.54e-02 
  2500| 4.94e-06  1.54e-04  8.21e-05  8.50e-01  8.50e-01  0.00e+00  1.60e-02 
  2600| 4.13e-06  1.29e-04  6.87e-05  8.50e-01  8.50e-01  0.00e+00  1.66e-02 
  2700| 3.46e-06  1.08e-04  5.75e-05  8.50e-01  8.50e-01  0.00e+00  1.72e-02 
  2760| 3.11e-06  9.69e-05  5.17e-05  8.50e-01  8.50e-01  0.00e+00  1.75e-02 
----------------------------------------------------------------------------
Status: Solved
Timing: Solve time: 1.76e-02s
    Lin-sys: nnz in L factor: 61, avg solve time: 3.53e-07s
    Cones: avg projection time: 5.11e-06s
----------------------------------------------------------------------------
Error metrics:
dist(s, K) = 2.2132e-09, dist(y, K*) = 9.3131e-10, s'y/m = -9.7499e-11
|Ax + s - b|_2 / (1 + |b|_2) = 3.1057e-06
|A'y + c|_2 / (1 + |c|_2) = 9.6854e-05
|c'x + b'y| / (1 + |c'x| + |b'y|) = 5.1706e-05
----------------------------------------------------------------------------
c'x = 0.8495, -b'y = 0.8497
============================================================================
0.8495352484237209

julia> @show MathProgBase.status(m)
MathProgBase.status(m) = :Optimal
:Optimal

julia> @show MathProgBase.getsolution(m)
MathProgBase.getsolution(m) = [0.2585352976879626,0.5909978259665712,0.8495352484237209,2.950007731763664,3.7708023318365167,4.2743511306364275e-8,3.27920973115259]
7-element Array{Float64,1}:
 0.258535  
 0.590998  
 0.849535  
 2.95001   
 3.7708    
 4.27435e-8
 3.27921   

julia> @show MathProgBase.getdual(m)
MathProgBase.getdual(m) = [-0.9999999998163523,0.08493642159298997,0.0,0.0,0.06370292305356898,0.0,0.06683902395033697,0.060756649730726565,-0.3656175766611347,0.02761388683088022,-0.23500428789197625,0.9999862689712925,0.0,0.0,0.0,0.027613898433779184,0.06944341154397124,-0.23500438167958973,0.34927248091608026,-0.8357844423855383,0.9999866469605179,0.0,0.0,0.0]
24-element Array{Float64,1}:
 -1.0      
  0.0849364
  0.0      
  0.0      
  0.0637029
  0.0      
  0.066839 
  0.0607566
 -0.365618 
  0.0276139
 -0.235004 
  0.999986 
  0.0      
  0.0      
  0.0      
  0.0276139
  0.0694434
 -0.235004 
  0.349272 
 -0.835784 
  0.999987 
  0.0      
  0.0      
  0.0    

[chriscoey]

Below is a second example where Mosek getdual fails. For this one, SCS getdual also fails (I will submit an issue to them):

c = [0.0,0.0,0.0,0.0,0.0,0.0,1.0,1.0]
b = [10.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,1.0,0.0,0.0,0.0,0.0]

I = [1,2,10,11,13,16,19,20,22,25,1,3,10,11,13,16,19,20,22,25,1,4,10,11,13,16,19,20,22,25,1,5,10,19,1,6,10,11,13,16,19,20,22,25,1,7,10,11,13,16,19,20,22,25,8,15,9,24]
J = [1,1,1,1,1,1,1,1,1,1,2,2,2,2,2,2,2,2,2,2,3,3,3,3,3,3,3,3,3,3,4,4,4,4,5,5,5,5,5,5,5,5,5,5,6,6,6,6,6,6,6,6,6,6,7,7,8,8]
V = [1.0,1.0,-0.44999999999999996,0.45000000000000007,-0.4500000000000001,5.551115123125783e-17,-0.44999999999999996,0.45000000000000007,-0.4500000000000001,5.551115123125783e-17,1.0,1.0,-0.7681980515339464,0.31819805153394654,-0.13180194846605373,5.551115123125783e-17,-0.7681980515339464,0.31819805153394654,-0.13180194846605373,5.551115123125783e-17,1.0,1.0,-0.9000000000000001,1.102182119232618e-16,-1.3497838043956718e-32,1.232595164407831e-32,-0.9000000000000001,1.102182119232618e-16,-1.3497838043956718e-32,1.232595164407831e-32,1.0,1.0,-0.22500000000000003,-0.22500000000000003,1.0,1.0,-0.11250000000000003,0.1125,-0.11249999999999999,-1.3877787807814457e-17,-0.11250000000000003,0.1125,-0.11249999999999999,-1.3877787807814457e-17,1.0,1.0,-8.436148777472949e-34,1.3777276490407724e-17,-0.22500000000000003,-1.5407439555097887e-33,-8.436148777472949e-34,1.3777276490407724e-17,-0.22500000000000003,-1.5407439555097887e-33,1.0,-1.0,1.0,-1.0]
A = sparse(I, J, V, length(b), length(c))

cone_con = [(:NonNeg,[1]),(:NonPos,[2,3,4,5,6,7,8,9]),(:SDP,10:15),(:Zero,16:18),(:SDP,19:24),(:Zero,25:27)]
cone_var = [(:NonNeg,[1,2,3,4,5,6,7,8])]

using MathProgBase, Mosek

m = MathProgBase.ConicModel(MosekSolver())
MathProgBase.loadproblem!(m, c, A, b, cone_con, cone_var)
MathProgBase.optimize!(m)

@show MathProgBase.status(m)
@show MathProgBase.getsolution(m)
@show MathProgBase.getdual(m)

[chriscoey]

another even smaller problem on which both Mosek and SCS fail with getdual:

c = [-0.0,-0.0,-0.0,-0.0,-0.0,-0.0,-1.0]
b = [10.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]

I = [1,2,8,9,10,11,1,3,8,9,10,11,1,4,8,9,10,11,1,5,8,1,6,8,9,10,11,1,7,8,9,10,11,8,10]
J = [1,1,1,1,1,1,2,2,2,2,2,2,3,3,3,3,3,3,4,4,4,5,5,5,5,5,5,6,6,6,6,6,6,7,7]
V = [1.0,1.0,-0.44999999999999996,0.45000000000000007,-0.4500000000000001,0.0,1.0,1.0,-0.7681980515339464,0.31819805153394654,-0.13180194846605373,0.0,1.0,1.0,-0.9000000000000001,0.0,0.0,0.0,1.0,1.0,-0.22500000000000003,1.0,1.0,-0.11250000000000003,0.1125,-0.11249999999999999,0.0,1.0,1.0,0.0,0.0,-0.22500000000000003,0.0,1.0,1.0]
A = sparse(I, J, V, length(b), length(c))

cone_con = [(:NonNeg,[1]),(:NonPos,[2,3,4,5,6,7]),(:SDP,8:10),(:Zero,11:11)]
cone_var = [(:NonNeg,[1,2,3,4,5,6]),(:Free,[7])]

using MathProgBase, Mosek

m = MathProgBase.ConicModel(MosekSolver())
MathProgBase.loadproblem!(m, c, A, b, cone_con, cone_var)
MathProgBase.optimize!(m)

@show MathProgBase.status(m)
@show MathProgBase.getsolution(m)
@show MathProgBase.getdual(m)

[mlubin]

@ulfworsoe, any idea what might be causing this crash when accessing duals for SDPs?

[ulfworsoe]

I can see you already figured it out. Works in the HEAD revision, so I'll close this.