In [2]:
using Convex
using SCS # I think I fixed it
#using ECOS
using Plots
using LinearAlgebra
In [3]:
include("plotregion.jl")
Out[3]:
Main.PlotRegion
In [4]:
A1 = [-2.0 1; 
      -1 2; 
       1 0]
b = [2.0; 7; 3]

A = [A1 Matrix{Float64}(I,3,3)] # form the problem with slacks.

PlotRegion.plotregion(A,b)
Out[4]:
In [5]:
# Convert the problem into standard form.
cs = [-1 -2 0 0 0]'
AS = A
Out[5]:
3×5 Matrix{Float64}:
 -2.0  1.0  1.0  0.0  0.0
 -1.0  2.0  0.0  1.0  0.0
  1.0  0.0  0.0  0.0  1.0
In [24]:
""" Solve the central-path problem for interior point methods. """
function ip_central(c,A,b,tau)
    x = Variable(length(c))
    constraints = Constraint[
        A*x == b
    ]
    p = minimize(c'*x - tau*sum(log(x)), constraints)
    
    solve!(p, SCS.Optimizer; silent=true)
    #solve!(p, ECOS.Optimizer; silent = true)
    return x.value, p
end
ip_central(cs,AS,b,10.0)[1]
Out[24]:
5×1 Matrix{Float64}:
 2.1944265108421326
 2.1971437784015126
 4.191709445576553
 4.800138904383412
 0.8055738457513615
In [7]:
ip_central(cs,AS,b,1e-7)[1]
Out[7]:
5×1 Matrix{Float64}:
 2.9999995053211297
 4.999999384661257
 2.9999996153139215
 7.2995288733118e-7
 4.927979362203644e-7
In [26]:
taus = vec([1000 100 10 7.5 5 3.5 2 1 0.75 0.5 0.35 0.20 10.0.^(range(-1,stop=-7,length=10))'])

p = PlotRegion.plotregion(AS,b)
for tau in taus
    x = ip_central(cs, AS, b, tau)[1]
    scatter!([x[1]],[x[2]],label="", marker_z=log10(tau))
end
p
Out[26]:
In [27]:
taus = vec([10 7.5 5 3.5 2 1 0.75 0.5 0.35 0.20 10.0.^(range(-1,stop=-7,length=10))'])

p = PlotRegion.plotregion(AS,b)
for tau in taus
    x = ip_central([-1 0 0 0 0.0]', AS, b, tau)[1]
    scatter!([x[1]],[x[2]],label="", marker_z=log10(tau))
end
p
Out[27]:
In [10]:
x = ip_central([-1,-2.0,0,0,0], AS, b, 0.0002)
@show x[1]

x[1] = [3.00000036770714; 4.999971118925543; 3.000029604328589; 5.83496815133341e-5; -2.567714595444165e-7;;]
Out[10]:
5×1 Matrix{Float64}:
  3.00000036770714
  4.999971118925543
  3.000029604328589
  5.83496815133341e-5
 -2.567714595444165e-7
In [11]:
A1 = [-2.0 1; 
      -1 2; 
       1 0]
b = [2.0; 7; 3]

AS = [A1 Matrix{Float64}(I,3,3)] # form the problem with slacks.
cs = [1 0 0 0 0]'
#cs = [-1 -2 0 0 0]'

taus = vec([10 7.5 5 3.5 2 1 0.75 0.5 0.35 0.20 10.0.^(range(-1,stop=-7,length=10))'])

p = PlotRegion.plotregion(AS,b)
for tau in taus
    x = ip_central(cs, AS, b, tau)[1]
    scatter!([x[1]],[x[2]],label="", color="red")
end
p
Out[11]:

Show a centered vs. uncentered step¶

In [12]:
# Convert the problem into standard form.
cs = [-1 -2 0 0 0]'
AS = A

tau = 1.0
x0,prob = ip_central(cs,AS,b,tau)
x = copy(x0)
x[1] += 0.5
x[2] += 0.5

lam = vec(prob.constraints[1].dual)
@show lam

# show the region and the starting point
plt = PlotRegion.plotregion(AS,b)
scatter!([x[1]],[x[2]],label="", color="red")

# compute the steps 
n = length(cs)
m = size(AS,1)

s = tau./x
J = [zeros(n,n) AS'         Matrix{Float64}(I,n,n); 
     AS        zeros(m,m)   zeros(m,n);
    Diagonal(vec(s))  zeros(n,m) Diagonal(vec(x))]

mu = dot(x,s)/n
sigma = 0.5

F = [s + AS'*lam - cs; AS*x - b; x.*s]

Fc = [s + AS'*lam .- cs; AS*x .- b; x.*s .- sigma*mu ]

p = J\-F
pc = J\-Fc

plot!([x[1];x[1] + p[1]], [x[2];x[2] + p[2]], label="Affine")
plot!([x[1];x[1] + pc[1]], [x[2];x[2] + pc[2]], label = "Centered")

plt

lam = [-0.3307127844906483, -0.9507406082818322, -2.9873654286696003]
Out[12]:
In [13]:
xf = [x; lam; s];

[xf+p xf+pc]
Out[13]:
13×2 Matrix{Float64}:
  2.89567      2.77989
  4.9671       4.63627
  2.82424      2.9235
 -0.0385367    0.507347
  0.10433      0.220114
 -0.0218177   -0.17632
 -0.985621    -0.967539
 -2.05616     -2.51661
  0.0269049    0.196429
 -0.00694125   0.111399
  0.0218177    0.17632
  0.985621     0.967539
  2.05616      2.51661
In [14]:
taus = vec([10 7.5 5 3.5 2 1 0.75 0.5 0.35 0.20 10.0.^(range(-1,stop=-7,length=10))'])

p = PlotRegion.plotregion(AS,b)
for tau in taus
    x = ip_central(cs, AS, b, tau)[1]
    @show tau
    @show x
    scatter!([x[1]],[x[2]],label="", color="red")
end
p

tau = 10.0
x = [2.1944265108421326; 2.1971437784015126; 4.191709445576553; 4.800138904383412; 0.8055738457513615;;]
tau = 7.5
x = [2.227159268446235; 2.3655281162143553; 4.088790491688571; 4.496103030490232; 0.7728407732588833;;]
tau = 5.0
x = [2.290049811835893; 2.696802563321335; 3.8832974702777743; 3.8964445445195715; 0.7099502781341624;;]
tau = 3.5
x = [2.3624744798306696; 3.075130982719054; 3.649818404167572; 3.212212267171172; 0.6375256619051859;;]
tau = 2.0
x = [2.497985999456187; 3.71257438082063; 3.283397706346886; 2.0728371924803395; 0.5020140874401607;;]
tau = 1.0
x = [2.66521989488419; 4.3067280411953295; 3.023711668119117; 1.0517638730932535; 0.33478004880715173;;]
tau = 0.75
x = [2.7259021239376233; 4.4715337854748745; 2.9802704196537193; 0.7828345806851613; 0.2740978407483049;;]
tau = 0.5
x = [2.798901652585561; 4.641529087482228; 2.9562742372819337; 0.515843465035907; 0.20109836432005826;;]
tau = 0.35
x = [2.850249210257786; 4.746122956682109; 2.9543754498205064; 0.35800330294821897; 0.1497507761304401;;]
tau = 0.2
x = [2.9087106579273687; 4.853024380665235; 2.9643969456418366; 0.20266189243507912; 0.09128935311475149;;]
tau = 0.1
x = [2.952217342042687; 4.925779971619681; 2.978654640282322; 0.1006574404484018; 0.047782593624758216;;]
tau = 0.021544346900318832
x = [2.9893424261822443; 4.983884214501555; 2.994800642422667; 0.02157400235405059; 0.010657586318152501;;]
tau = 0.004641588833612777
x = [2.997461609037214; 4.996429047180418; 2.9984937609149123; 0.004603711429527302; 0.0025379834484078314;;]
tau = 0.001
x = [2.9991946161723657; 4.999151743432025; 2.9992368020410116; 0.0008915922472728816; 0.0008048294321981193;;]
tau = 0.00021544346900318845
x = [2.9997465473938543; 4.999690348496971; 2.9998027461468317; 0.0003656559608840747; 0.0002532724116719394;;]
tau = 4.641588833612782e-5
x = [3.000019789752743; 5.00003099451976; 3.00000849664363; -4.2018861696248565e-5; -1.974908477987863e-5;;]
tau = 1.0e-5
x = [3.0000125267915587; 4.9999972744355645; 3.000027647697383; 1.7953616662693114e-5; -1.2502224304444305e-5;;]
tau = 2.1544346900318822e-6
x = [3.000000937409824; 5.000000088052242; 3.0000017767205414; 7.686406690782814e-7; -9.338808988456467e-7;;]
tau = 4.641588833612782e-7
x = [2.9999997157547327; 4.999999563783239; 2.999999857511467; 5.851101773751007e-7; 2.8328491662582974e-7;;]
tau = 1.0e-7
x = [2.9999995053211297; 4.999999384661257; 2.9999996153139215; 7.2995288733118e-7; 4.927979362203644e-7;;]
Out[14]:
In [15]:
x = ip_central([-1,-2.0,0,0,0], AS, b, 0.0002)
@show x[1]

x[1] = [3.00000036770714; 4.999971118925543; 3.000029604328589; 5.83496815133341e-5; -2.567714595444165e-7;;]
Out[15]:
5×1 Matrix{Float64}:
  3.00000036770714
  4.999971118925543
  3.000029604328589
  5.83496815133341e-5
 -2.567714595444165e-7