Can I have a faster code?

[Giacomo Canciello]

Hello, my problem is:
min norm(AXB-C)
subjact to
X>=0
where X is a diagonal block variable
After 100 iteration the result is excellent, but the running time of the program is 7 hours.
Can I have a faster code?thanks

function RR=minnew(A,B,C)

R1=sdpvar(2,2);
R2=sdpvar(2,2);
R3=sdpvar(2,2);
R4=sdpvar(2,2);
R5=sdpvar(2,2);
R6=sdpvar(2,2);
R7=sdpvar(2,2);
R8=sdpvar(2,2);
R9=sdpvar(2,2);
R10=sdpvar(2,2);
R11=sdpvar(2,2);
R12=sdpvar(2,2);
X=blkdiag(R1,R2,R3,R4,R5,R6,R7,R8,R9,R10,R11,R12);
obj=norm(A*X*B-C);
v=[X>=0];
opt=sdpsettings('solver','sedumi','sedumi.eps',1e-20,'sedumi.maxiter',100);
solvesdp(v,obj,opt);
RR=double(X);

[Erling Andersen]

Have tried MOSEK (or another optimizer) instead of SeDuMi.

Since I do not have the data I cannot do the expiriment myself.

[Giacomo Canciello]

If I use mosek the result is NaN.Why?

obj=norm(A*X*B-C);
v=[X>=0];
opt=sdpsettings('solver','mosek');
solvesdp(v,obj,opt);
RR=double(X);

[Erling Andersen]

That must be caused by a bug in Yalmip or MOSEK. 

At

   http://users.isy.liu.se/johanl/yalmip/pmwiki.php?n=Solvers.MOSEK

it is described how you can dump the problem to disk that you can email to sup...@mosek.com.

Then we shall investigate the issue.

 

[Johan Löfberg]

What are the dimensions of A, B and C? Real or complex?

[Johan Löfberg]

Do you have mosek installed?

What diagnostic code is returned by solvesdp (the field 'problem')

Try with debug turned on

opt=sdpsettings('solver','mosek','debug',1);
solvesdp(v,obj,opt);


[Giacomo Canciello]

Sorry ,now Mosek is installed.Sorry,I thought that matlab advise me.
Now I 'm trying mosek,I'll know when it ends.

However the dimension is:
A:34996x24 real
B:24x2 real
C:34996x2 real

P.S. If i use LMIlab of robust control the running time is only 2-3 min,but the results are worse that yalmip because i can't declare a constraint semidefinite X>=0 but only definite X>0.

[Johan Löfberg]

There is no difference between >0 and >= in practice. Most solvers use finite-precision algorithms which approach the solution from possibly infeasible directions, and you can even get slightly infeasible solutions. If you want to have a strict solution, you constrain it to >= eps*I where eps is a constant you pick, which is sensible in terms of solver precision, machine precision, and application background. Additionally, saying >0 is problematic since there is no minimizer (only an infimizer)

You should probably change the measure you minimize. The norm (2-norm here) requires an SDP of size 34998 x 34998 (I am surprised the code even runs, and very very surprised when you say lmilab solves the problem in 2 minutes. Care to share the code?)

[Giacomo Canciello]

I use all data in LMIlab (A:34996x24 B:24x2 C:34996x2), on the contrary I use 1 sample every 12 (f(1:12:end)) (A:2920x24 B:24x2 C:2920x2) because it's block.

In LMIlab I turned my problem in the following way:
min(AXB-C) is equal to (AXB-C)' (AXB-C)=epsilon*I -->(B'XA')(AXB-C)=epsilon*I
maxeig(B'XA'AXB-B'XA'C-C'AXB+C'C-epsion^2*I)>0
B'XA'AXB-B'XA'C-C'AXB+C'C-epsion^2*I<0
This can be written as complment schur:
[C'*C'-C'*A*X*B-B'*X*A'*C-epsilon^2*I      -B'*X;    -X*B     -inv(A'*A) ] < 0

[C'*C'-C'*A*X*B-B'*X*A'*C     -B'*X;    -X*B     -inv(A'*A) ] < epsilon^2 [I 0;0 0]
finally i solve with gevp

function [Xopt,errore]=minLMI(A,B,C)

[m,r]=size(C);
[m,n]=size(A);
nm = n/r;
ee=10^-9;

setlmis([])
X=lmivar(1,repmat([r,1],nm,1));
Y=lmivar(1,[r,1]);


lmiterm([-1,1,1,X],1,1);
lmiterm([1,1,1,0],-ee*eye(nm*r))

%A<[Y 0;0 0]
lmiterm([2 1 1 X],-C'*A,B,'s')
lmiterm([2,1,1,0],C'*C)
lmiterm([2,1,2,X],-B',1)
lmiterm([2,2,2,0],-inv(A'*A))
lmiterm([-2,1,1,Y],1,1)
lmiterm([-2,2,2,0],zeros(n,r))

%0<B1
lmiterm([-3,1,1,0],eye(r))

%Y<lambda*B1
lmiterm([4,1,1,Y],1,1)
lmiterm([-4,1,1,0],eye(r))


LMIsist=getlmis;

options=([1e-15,500,1e10,30,0]);

[errmin,xvalutato]=gevp(LMIsist,1,options);
Xopt=dec2mat(LMIsist,xvalutato,X);

errore=sqrt(errmin)

P.S. If i use thi structure with yalmip tre result is horrible:

function RR=minnew3(A,B,C)
[m,r]=size(C);
[m,n]=size(A);
nm = n/r;
ee=sdpvar(1,1);

R1=sdpvar(2,2);
R2=sdpvar(2,2);
R3=sdpvar(2,2);
R4=sdpvar(2,2);
R5=sdpvar(2,2);
R6=sdpvar(2,2);
R7=sdpvar(2,2);
R8=sdpvar(2,2);
R9=sdpvar(2,2);
R10=sdpvar(2,2);
R11=sdpvar(2,2);
R12=sdpvar(2,2);
X=blkdiag(R1,R2,R3,R4,R5,R6,R7,R8,R9,R10,R11,R12);

F=[[C'*C-C'*A*X*B-B'*X*A'*C-ee*eye(r)   -B'*X;
               -X*B                -inv(A'*A) ] < 0,X>=0];

sol=solvesdp(F,ee);
RR=double(X);

[Giacomo Canciello]

sorry, there is a mistake in my post:maxeig(B'XA'AXB-B'XA'C-C'AXB+C'C-epsion^2*I)<0

[Johan Löfberg]

Why are you solving a GEVP? It is a simple linear SDP

YALMIP performs Schur complements etc too to represent the norm, but does not know about the low-rank structure which you exploit when you pull out A'*A. Hence, you should model the norm manually in this case to get a small SDP.

To proceed, you would have to post the data.

[Giacomo Canciello]

I must change the norm problem with LMI problem with semidefinite bond if i want to use LMIlab, unlike I use norm problem in yalmip.
Do you want A,B, C data?

[Johan Löfberg]

Yes.

[Giacomo Canciello]

I am attaching both the data complete (datacomplete) and data 1 sample every 12 (data12)

[Johan Löfberg]

Your model is a numerical catastrophe. The inverse of A'A makes no sense since the matrix is singular

>> eig(A'*A)

ans =

     1.468926470288500e-15

     1.468926470310432e-15

     4.166500229953062e-13

     4.166500229953183e-13

     4.844656278979242e-13

     4.844656278979328e-13

     5.797188272585414e-13

     5.797188272585999e-13

     9.808581270744462e-13

     9.808581270744545e-13

     1.702620885435675e-12

     1.702620885435733e-12

     2.833219150567907e-12

     2.833219150567916e-12

     4.372633994054583e-12

     4.372633994054755e-12

     1.206465435729847e-11

     1.206465435729873e-11

     2.450315844979118e-11

     2.450315844979123e-11

     1.657671493142084e-10

     1.657671493142087e-10

     1.041004009811770e-09

     1.041004009811770e-09

[Johan Löfberg]

...and the data is basically numerical noise. C has all elements on the order of 10^-7, and A 10^-12. You have to clean up you model to be able to get something of value from this. By brutally scaling the data to somewhat more reasonable numbers, and ad-hoc regularizing your inverse, i.e. replacing A'A with A'A+1e-5*eye(24), you can get better solutions than the value lmilab happens to spit out when encountering this numerical disaster 

RR = minnew3(A*1e5,B,C*1e5);norm(A*RR*B-C)

ans =

     4.716288311534027e-05

[Giacomo Canciello]

very very thanks!!!