Biostat M280: Compressed Sensing by Linear Programming (LP)

Adapted from candes.m by Dr. Hua Zhou, Mar 3, 2016

Display system information.

In [1]:
versioninfo()
Julia Version 0.4.3
Commit a2f713d (2016-01-12 21:37 UTC)
Platform Info:
  System: Darwin (x86_64-apple-darwin13.4.0)
  CPU: Intel(R) Core(TM) i7-3720QM CPU @ 2.60GHz
  WORD_SIZE: 64
  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Sandybridge)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.3

Generate a sparse signal and sub-sampling

In [2]:
using Gadfly

# random seed
srand(20160303)
# Size of signal
n = 1024
# Sparsity (# nonzeros) in the signal
s = 20
# Number of samples (undersample by a factor of 8) 
m = 128

# Generate and display the signal
x0 = zeros(n)
x0[rand(1:n, s)] = randn(s)
# Generate the random sampling matrix
A = randn(m, n) / m
# Subsample by multiplexing
y = A * x0

# plot the true signal
plot(x=1:n, y=x0, Geom.line, Guide.title("True signal x_0"))
Out[2]:
x -2000 -1500 -1000 -500 0 500 1000 1500 2000 2500 3000 3500 -1500 -1450 -1400 -1350 -1300 -1250 -1200 -1150 -1100 -1050 -1000 -950 -900 -850 -800 -750 -700 -650 -600 -550 -500 -450 -400 -350 -300 -250 -200 -150 -100 -50 0 50 100 150 200 250 300 350 400 450 500 550 600 650 700 750 800 850 900 950 1000 1050 1100 1150 1200 1250 1300 1350 1400 1450 1500 1550 1600 1650 1700 1750 1800 1850 1900 1950 2000 2050 2100 2150 2200 2250 2300 2350 2400 2450 2500 2550 2600 2650 2700 2750 2800 2850 2900 2950 3000 -2000 0 2000 4000 -1500 -1400 -1300 -1200 -1100 -1000 -900 -800 -700 -600 -500 -400 -300 -200 -100 0 100 200 300 400 500 600 700 800 900 1000 1100 1200 1300 1400 1500 1600 1700 1800 1900 2000 2100 2200 2300 2400 2500 2600 2700 2800 2900 3000 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 -8.0 -7.8 -7.6 -7.4 -7.2 -7.0 -6.8 -6.6 -6.4 -6.2 -6.0 -5.8 -5.6 -5.4 -5.2 -5.0 -4.8 -4.6 -4.4 -4.2 -4.0 -3.8 -3.6 -3.4 -3.2 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0 5.2 5.4 5.6 5.8 6.0 6.2 6.4 6.6 6.8 7.0 -10 -5 0 5 10 -8.0 -7.5 -7.0 -6.5 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 y True signal x_0

Solve LP by calling Gurobi directly (Gurobi.jl)

Gurobi model formulation:
objective: (1/2) x' H x + f' x
s.t. A x <= b
Aeq x = beq
lb <= x <= ub
Refer to Gurobi.jl documentation for setting up the model.

In [3]:
using Gurobi
env = Gurobi.Env()
setparams!(env, OutputFlag=1) # display log

# Construct the model
model = gurobi_model(env;
    name = "cs",
    f = ones(2 * n),
    Aeq = [A -A],
    beq = y,
    lb = zeros(2 * n))

# Run optimization
optimize(model)

# Show results
sol = get_solution(model)
xsol = sol[1:n] - sol[n + 1:end]

plot(x=1:n, y=x0, Geom.point)
plot(x=1:n, y=xsol, Geom.line, Guide.title("Reconstructed signal overlayed with x0"))
Optimize a model with 128 rows, 2048 columns and 262144 nonzeros
Coefficient statistics:
  Matrix range    [7e-08, 3e-02]
  Objective range [1e+00, 1e+00]
  Bounds range    [0e+00, 0e+00]
  RHS range       [6e-04, 1e-01]

Concurrent LP optimizer: dual simplex and barrier
Showing barrier log only...

Presolve time: 0.09s
Presolved: 128 rows, 2048 columns, 262144 nonzeros

Ordering time: 0.00s

Barrier statistics:
 AA' NZ     : 8.128e+03
 Factor NZ  : 8.256e+03 (roughly 1 MByte of memory)
 Factor Ops : 7.073e+05 (less than 1 second per iteration)
 Threads    : 3

                  Objective                Residual
Iter       Primal          Dual         Primal    Dual     Compl     Time
   0   1.00460219e+03  0.00000000e+00  3.13e-14 0.00e+00  7.20e-01     0s
   1   2.56436546e+02  6.71121241e+00  4.26e-14 4.44e-16  1.22e-01     0s
   2   4.47033257e+01  1.15222875e+01  1.14e-14 4.44e-16  1.62e-02     0s
   3   2.35369766e+01  1.56681328e+01  4.83e-13 4.44e-16  3.84e-03     0s
   4   1.97642318e+01  1.71461269e+01  3.43e-13 4.44e-16  1.28e-03     0s
   5   1.89622258e+01  1.76971372e+01  3.11e-13 8.88e-16  6.18e-04     0s
   6   1.84876234e+01  1.81452336e+01  1.88e-12 4.44e-16  1.67e-04     0s
   7   1.83880675e+01  1.81830002e+01  1.81e-12 4.44e-16  1.00e-04     0s
   8   1.82351093e+01  1.81981958e+01  2.37e-13 1.29e-06  1.80e-05     0s
   9   1.82247106e+01  1.82245859e+01  3.45e-12 1.22e-15  6.09e-08     0s
  10   1.82246148e+01  1.82246147e+01  9.56e-13 1.89e-15  6.09e-11     0s

Barrier solved model in 10 iterations and 0.20 seconds
Optimal objective 1.82246148e+01

Crossover log...

     108 DPushes remaining with DInf 0.0000000e+00                 0s

Solved with dual simplex
Solved in 308 iterations and 0.25 seconds
Optimal objective  1.822461471e+01
Out[3]:
x -2000 -1500 -1000 -500 0 500 1000 1500 2000 2500 3000 3500 -1500 -1450 -1400 -1350 -1300 -1250 -1200 -1150 -1100 -1050 -1000 -950 -900 -850 -800 -750 -700 -650 -600 -550 -500 -450 -400 -350 -300 -250 -200 -150 -100 -50 0 50 100 150 200 250 300 350 400 450 500 550 600 650 700 750 800 850 900 950 1000 1050 1100 1150 1200 1250 1300 1350 1400 1450 1500 1550 1600 1650 1700 1750 1800 1850 1900 1950 2000 2050 2100 2150 2200 2250 2300 2350 2400 2450 2500 2550 2600 2650 2700 2750 2800 2850 2900 2950 3000 -2000 0 2000 4000 -1500 -1400 -1300 -1200 -1100 -1000 -900 -800 -700 -600 -500 -400 -300 -200 -100 0 100 200 300 400 500 600 700 800 900 1000 1100 1200 1300 1400 1500 1600 1700 1800 1900 2000 2100 2200 2300 2400 2500 2600 2700 2800 2900 3000 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 -8.0 -7.8 -7.6 -7.4 -7.2 -7.0 -6.8 -6.6 -6.4 -6.2 -6.0 -5.8 -5.6 -5.4 -5.2 -5.0 -4.8 -4.6 -4.4 -4.2 -4.0 -3.8 -3.6 -3.4 -3.2 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0 5.2 5.4 5.6 5.8 6.0 6.2 6.4 6.6 6.8 7.0 -10 -5 0 5 10 -8.0 -7.5 -7.0 -6.5 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 y Reconstructed signal overlayed with x0

Solve LP by DCP (disciplined convex programming) interface Convex.jl

In [5]:
using Convex

# Use Gurobi solver
using Gurobi
solver = GurobiSolver(OutputFlag=1)
set_default_solver(solver)

## Use Mosek solver
#using Mosek
#solver = MosekSolver(LOG=1)
#set_default_solver(solver)

## Use SCS solver
#using SCS
#solver = SCSSolver(verbose=1)
#set_default_solver(solver)

# Set up optimizaiton problem
x = Variable(n)
problem = minimize(norm(x, 1))
problem.constraints += A * x == y

# Solve the problem
@time solve!(problem)

# Display the solution
plot(x=1:n, y=x0, Geom.point)
plot(x=1:n, y=xsol, Geom.line, Guide.title("Reconstructed signal overlayed with x0"))
Optimize a model with 2177 rows, 2049 columns and 136193 nonzeros
Coefficient statistics:
  Matrix range    [7e-08, 1e+00]
  Objective range [1e+00, 1e+00]
  Bounds range    [0e+00, 0e+00]
  RHS range       [6e-04, 1e-01]

Concurrent LP optimizer: dual simplex and barrier
Showing barrier log only...

Presolve removed 1025 rows and 1 columns
Presolve time: 0.06s
Presolved: 1152 rows, 2048 columns, 133120 nonzeros

Ordering time: 0.00s

Barrier statistics:
 Free vars  : 1024
 AA' NZ     : 1.392e+05
 Factor NZ  : 1.590e+05 (roughly 3 MBytes of memory)
 Factor Ops : 2.343e+07 (less than 1 second per iteration)
 Threads    : 3

                  Objective                Residual
Iter       Primal          Dual         Primal    Dual     Compl     Time
   0   1.01141109e+03  0.00000000e+00  7.11e-15 1.00e+00  6.17e-01     0s
   1   4.03536450e+02  6.25798792e+00  6.29e-08 1.69e-03  1.94e-01     0s
   2   1.11879208e+02  1.08616096e+01  1.70e-08 4.64e-04  4.93e-02     0s
   3   4.10024083e+01  1.49413770e+01  1.94e-08 1.35e-04  1.27e-02     0s
   4   2.55415844e+01  1.68876683e+01  1.26e-08 6.62e-05  4.23e-03     0s
   5   2.12370478e+01  1.76323916e+01  1.17e-08 3.63e-05  1.76e-03     0s
   6   1.92591230e+01  1.78986854e+01  6.90e-09 1.65e-05  6.64e-04     0s
   7   1.85725233e+01  1.81726688e+01  6.25e-09 6.03e-06  1.95e-04     0s
   8   1.84775909e+01  1.81487061e+01  4.51e-09 5.60e-06  1.61e-04     0s
   9   1.82954753e+01  1.82136656e+01  2.36e-09 1.17e-06  4.00e-05     0s
  10   1.82252231e+01  1.82237444e+01  2.15e-10 6.95e-08  7.22e-07     0s
  11   1.82246153e+01  1.82246138e+01  3.31e-10 3.15e-11  7.22e-10     0s
  12   1.82246147e+01  1.82246147e+01  7.34e-11 1.17e-14  1.13e-12     0s

Barrier solved model in 12 iterations and 0.21 seconds
Optimal objective 1.82246147e+01

Crossover log...

       0 DPushes remaining with DInf 6.7458121e-13                 0s

       0 PPushes remaining with PInf 0.0000000e+00                 0s

  Push phase complete: Pinf 0.0000000e+00, Dinf 6.7458121e-13      0s


Solved with dual simplex
Solved in 1319 iterations and 0.28 seconds
Optimal objective  1.822461471e+01
  
Out[5]:
x -2000 -1500 -1000 -500 0 500 1000 1500 2000 2500 3000 3500 -1500 -1450 -1400 -1350 -1300 -1250 -1200 -1150 -1100 -1050 -1000 -950 -900 -850 -800 -750 -700 -650 -600 -550 -500 -450 -400 -350 -300 -250 -200 -150 -100 -50 0 50 100 150 200 250 300 350 400 450 500 550 600 650 700 750 800 850 900 950 1000 1050 1100 1150 1200 1250 1300 1350 1400 1450 1500 1550 1600 1650 1700 1750 1800 1850 1900 1950 2000 2050 2100 2150 2200 2250 2300 2350 2400 2450 2500 2550 2600 2650 2700 2750 2800 2850 2900 2950 3000 -2000 0 2000 4000 -1500 -1400 -1300 -1200 -1100 -1000 -900 -800 -700 -600 -500 -400 -300 -200 -100 0 100 200 300 400 500 600 700 800 900 1000 1100 1200 1300 1400 1500 1600 1700 1800 1900 2000 2100 2200 2300 2400 2500 2600 2700 2800 2900 3000 -9 -8 -7 -6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 -8.0 -7.8 -7.6 -7.4 -7.2 -7.0 -6.8 -6.6 -6.4 -6.2 -6.0 -5.8 -5.6 -5.4 -5.2 -5.0 -4.8 -4.6 -4.4 -4.2 -4.0 -3.8 -3.6 -3.4 -3.2 -3.0 -2.8 -2.6 -2.4 -2.2 -2.0 -1.8 -1.6 -1.4 -1.2 -1.0 -0.8 -0.6 -0.4 -0.2 0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 2.4 2.6 2.8 3.0 3.2 3.4 3.6 3.8 4.0 4.2 4.4 4.6 4.8 5.0 5.2 5.4 5.6 5.8 6.0 6.2 6.4 6.6 6.8 7.0 -10 -5 0 5 10 -8.0 -7.5 -7.0 -6.5 -6.0 -5.5 -5.0 -4.5 -4.0 -3.5 -3.0 -2.5 -2.0 -1.5 -1.0 -0.5 0.0 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 5.5 6.0 6.5 7.0 y Reconstructed signal overlayed with x0
0.312172 seconds (8.56 k allocations: 21.860 MB, 1.17% gc time)
In [ ]: