Solver status is Infeasible


(Prateek Nallabolu) #1

Hello,
I am trying to solve a basis pursuit algorithm (Compressed sensing) using CVX. The algorithm is given as:
image
In my case, “A” is the multiplication of two matrices “Directivity” and “Baseband”. Since each element of the matrix “baseband” is also a complex 1-D array (1x1000), I couldn’t set up the algorithm directly in CVX. I used the “expression” feature in CVX to manually calculate the product “Ax”. The matrix “y” is of the size 10x1000.

cvx_begin 
variable x(399)
expression A(10,1000)
for ii=1:10
    for jj=1:19
        for kk=1:21
            A(ii,:)= A(ii,:)+ x((jj-1)*21+kk).*directivity(ii,kk).*baseband(jj,:);
        end
    end
end
minimize(norm(x,1))
subject to
y == A
cvx_end

Calling SDPT3 4.0: 20798 variables, 798 equality constraints
For improved efficiency, SDPT3 is solving the dual problem.

num. of constraints = 798
dim. of socp var = 798, num. of socp blk = 399
dim. of free var = 20000 *** convert ublk to lblk


SDPT3: Infeasible path-following algorithms


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

0|0.000|0.000|3.9e-01|4.2e+02|3.2e+09|-5.509656e-07 0.000000e+00| 0:0:20| spchol 1 1
1|1.000|0.984|2.6e-06|7.1e+00|1.1e+07|-3.422388e+04 -9.377909e+01| 0:0:38| spchol 1 1
2|0.082|0.683|2.6e-06|2.4e+00|8.2e+05|-7.236874e+04 -1.402673e+04| 0:0:55| spchol 1 1
3|0.326|0.157|3.8e-06|2.7e+00|6.7e+05|-5.367168e+05 -1.289798e+04| 0:1:12| spchol 1 1
4|1.000|0.289|3.3e-06|1.9e+00|8.1e+05|-3.782611e+06 -1.056407e+04| 0:1:26| spchol 1 1
5|1.000|0.027|3.2e-06|2.5e+00|2.4e+08|-5.820550e+08 -1.049181e+04| 0:1:41| spchol 1 1
6|0.083|0.254|5.0e-06|1.9e+00|5.5e+07|-1.133612e+09 -2.776567e+05| 0:1:56| spchol 2 2
7|0.159|0.026|2.6e-05|2.4e+00|6.1e+10|-1.735589e+11 -2.842752e+05| 0:2:10| spchol
linsysolve: Schur complement matrix not positive definite
switch to LU factor. splu 22 30
8|0.003|0.001|5.1e-03|3.1e+00|1.3e+11|-1.804904e+11 -5.542437e+05| 0:2:26| splu 20 12
9|0.004|0.003|5.9e-02|3.7e+00|1.9e+11|-1.837358e+11 -1.394131e+06| 0:2:40| splu 17 ^18
10|0.001|0.001|3.9e-01|4.3e+00|2.6e+11|-1.840341e+11 -2.264338e+06| 0:2:55| splu 17 9
11|0.000|0.000|4.8e-01|4.9e+00|3.2e+11|-1.841321e+11 -2.666387e+06| 0:3:09| splu 30 ^16
12|0.000|0.000|6.2e-01|5.5e+00|3.8e+11|-1.841404e+11 -2.757989e+06| 0:3:24| splu 12 ^ 7
13|0.000|0.000|7.1e-01|6.2e+00|4.5e+11|-1.841485e+11 -2.827394e+06| 0:3:39| splu 11 23
14|0.000|0.000|7.2e-01|6.8e+00|5.1e+11|-1.841538e+11 -2.911517e+06| 0:3:53| splu 30 30
15|0.000|0.000|7.2e-01|7.4e+00|5.7e+11|-1.841541e+11 -2.934834e+06| 0:4:08| splu 12 ^17
16|0.000|0.000|7.4e-01|8.1e+00|6.4e+11|-1.841556e+11 -2.998431e+06| 0:4:23| splu 20 ^15
17|0.000|0.000|7.5e-01|8.7e+00|7.0e+11|-1.841571e+11 -3.066418e+06| 0:4:37| splu 25 ^14
18|0.000|0.000|7.7e-01|9.3e+00|7.7e+11|-1.841571e+11 -3.127167e+06| 0:4:53| splu 18 ^11
19|0.000|0.000|7.7e-01|9.9e+00|8.3e+11|-1.841574e+11 -3.120555e+06| 0:5:07| splu 16 12
20|0.000|0.000|7.8e-01|1.1e+01|8.9e+11|-1.841593e+11 -3.320728e+06| 0:5:22| splu 30 ^ 6
21|0.000|0.000|8.0e-01|1.1e+01|9.6e+11|-1.841600e+11 -3.522938e+06| 0:5:36| splu 12 ^21
22|0.000|0.000|8.1e-01|1.2e+01|1.0e+12|-1.841606e+11 -3.706480e+06| 0:5:51| splu 20 18
23|0.000|0.000|8.2e-01|1.2e+01|1.1e+12|-1.841612e+11 -3.901598e+06| 0:6:06| splu 28 ^12
24|0.000|0.000|8.3e-01|1.3e+01|1.1e+12|-1.841613e+11 -3.947350e+06| 0:6:21| splu 13 ^18
25|0.000|0.000|8.4e-01|1.4e+01|1.2e+12|-1.841616e+11 -3.976903e+06| 0:6:35| splu 18 ^23
26|0.000|0.000|8.6e-01|1.4e+01|1.3e+12|-1.841620e+11 -3.990878e+06| 0:6:50| splu 24 16
27|0.000|0.000|8.6e-01|1.5e+01|1.3e+12|-1.841627e+11 -4.328514e+06| 0:7:05| splu 30 11
28|0.000|0.000|8.8e-01|1.6e+01|1.4e+12|-1.841634e+11 -4.542460e+06| 0:7:19| splu 12 ^ 6
29|0.000|0.000|9.0e-01|1.6e+01|1.5e+12|-1.841640e+11 -4.642593e+06| 0:7:35| splu 13 ^ 7
30|0.000|0.000|9.1e-01|1.7e+01|1.5e+12|-1.841648e+11 -4.772688e+06| 0:7:49| splu 15 14
31|0.000|0.000|9.1e-01|1.7e+01|1.6e+12|-1.841650e+11 -4.917647e+06| 0:8:04|
sqlp stop: dual problem is suspected of being infeasible

number of iterations = 31
residual of dual infeasibility
certificate X = 1.50e-10
reldist to infeas. <= 2.82e-12
Total CPU time (secs) = 483.76
CPU time per iteration = 15.61
termination code = 2
DIMACS: 9.5e+00 0.0e+00 9.0e+02 0.0e+00 -1.0e+00 8.7e+00


Status: Infeasible
Optimal value (cvx_optval): +Inf

The status of the solver is Infeasible. I am not sure if I am using the “expression” feature correctly. I went through the CVX user guide but I couldn’t figure out the error. Can anyone please point out the error in the code?
I am not very familiar with CVX and any help will be greatly appreciated.


(Mark L. Stone) #2

I don’t understand what you have done, and can’t say whether it is correct, but it looks rather dubious to me. Why is x 399 by 1 rather than 1000 by 1?

Perhaps A could be a 3-D array, 10 by 1000 by 1000. Then you could specify the following as the constraints in CVX?

for i=1:1000
  y(:,i) == A(:,:,i)*x
end

Does that specify what you want, presuming you construct A correctly?

Anyhow, i don’t know much about compressed sensing, but the constraint(s) y == Ax (for instance per my for loop above) seems like it would not be the “correct” constraint to have, and would indeed be overconstrained such that it would be almost impossible to be feasible. Do you perhaps want a norm constraint on y-Ax or as a term in the objective function?


(Prateek Nallabolu) #3

Thank you for the suggestion of using a for loop in the constraints and making “A” 3-D array. I’ve restructured the CVX program without using expressions and used a for loop in the constraints. I am able to solve the algorithm now.

A = zeros(1000,399,10);
for ii=1:10
for jj=1:1000
for kk=1:399
col = uint8(mod(kk-1,21));
col1 = uint8(floor((kk-1)/21));
A(jj,kk,ii) = directivity(ii,col+1)*baseband(col1+1,jj);
end
end
end
cvx_begin
variable x(399,1)
minimize(norm(x,1))
subject to
for ii=1:10
for jj=1:1000
y(ii,jj) == A(jj,:,ii)*x
end
end
cvx_end

Sorry for not being clear with my problem earlier. I’ve tried to explain my constraint in the image below. Hope it explains why x is 399 by 1 and not 1000 by 1.

The equality constraint y=Ax works because I am considering a noiseless environment. In the case of noise, as you mentioned I will have a constraint on the L2 norm of y-Ax.
Thank you once again for your help.


(Mark L. Stone) #4

For the record, I corrected y(i,:) to y(:,i) in my previous post.


(Prateek Nallabolu) #5

The first approach using “expression” also works. I made a small mistake in the MATLAB program.
Thank you once again for your support.